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Abstract 

We consider disk accretion onto a magnetic star under the assump- 
tion that the stehar magnetic field permeates the disk and that the 
^ magnetosphere that hes between the disk and the star is force free. 

■ Using simphfied axisymmetric models (both semianalytic and numer- 

c/3 . ical), we study the time evolution of the magnetic field configuration 

^ ' induced by the relative rotation between the disk and the star. We 

show that, if both the star and the magnetosphere can be approxi- 
^ . mated as being perfectly conducting, then the behavior of the twisted 

H ' field lines depends on the magnitude of the surface conductivity of 

the disk. For any given relative azimuthal speed /S.v^ between the 
disk and the star (measured at the disk surface), there is a maximum 
surface conductivity Smax ~ such that, if the actual surface 

conductivity is smaller than Umax, then a steady-state field configu- 
ration can be established, whereas for larger values no steady state is 
possible and the field lines inflate and open up when a critical twist 
angle (which for an initially dipolar field is ~ 2 rad) is attained. We 
argue that thin astrophysical disks are likely to have surface conduc- 
tivities that exceed the local Smax except in regions where \5vA is 
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particularly small (such as the immediate vicinity of the corotation 
radius in a Keplerian disk). 

We find that, if the disk conductivity is high enough for the twist 
angle to continue growing until the field lines open, then, except un- 
der rather special circumstances, the radial magnetic field at the disk 
surface will grow to a level for which the radial migration of the field 
lines in the disk cannot be ignored. A similar conclusion was reached 
by Bardou & Heyvaerts (1996). We demonstrate, however, that the 
radial diffusion in the disk is much slower than the field-line expansion 
in the magnetosphere, which suggests that, contrary to the claim by 
Bardou & Heyvaerts, the twisting process does not result in a rapid 
expulsion of the field lines from the disk. 

We calculate the magnetospheric density redistribution brought 
about by the field-line expansion and show that inertial effects in the 
magnetosphere will slow the expansion (and invalidate the assumption 
of a force-free field) before the critical twist angle is reached. We ar- 
gue, however, that these effects merely slow down the opening process 
but do not inhibit it altogether. We find that the field opening dras- 
tically reduces the density near the apex of the expanding field lines 
(which typically elongate in a direction making an angle > 60° to the 
rotation axis) while also creating a pronounced density enhancement 
near the rotation axis. The former effect is conducive to the trigger- 
ing of microinstabilities (e.g. , ion-acoustic) in the current layer that 
forms along the direction of elongation, whereas the latter is evidently 
related to a mechanism for the formation of an axially outflowing con- 
densation that was previously identified in axisymmetric numerical 
simulations of such systems. 

We examine the possibility that the expanding field lines reconnect 
across the current layer before they open up, so that the twisting leads 
to a quasi-periodic process of field-line expansion and reconnection and 
not just to a one-time opening event. We tentatively conclude that 
hyperresistivity associated with tearing-mode turbulence could lead to 
efficient reconnection, but we speculate that even faster reconnection 
might be brought about by 3-D kinking of the twisted field lines that 
gave rise to thin current sheets. This question could be further pursued 
through MHD stability calculations and 3-D numerical simulations. 

Subject headings: accretion, accretion disks — MHD — stars: formation — 
stars: magnetic fields — stars: winds, outflows — stars: pre-main-sequence 
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1 Introduction 



The process of disk accretion onto a magnetic star is thought to have impor- 
tant consequences for the spin cvohition and observational characteristics of 
neutron stars, white dwarfs, and young stellar objects (YSOs). In particular, 
if the stellar magnetic field lines penetrate the disk, then they may transmit 
torques between the disk and the star. Such an interaction has been invoked 
to explain the spin-up and spin-down episodes observed in X-ray pulsars (e.g., 
Ghosh & Lamb 1978, 1979a, 1979b, hereafter GL; Wang 1987; Lovelace, Ro- 
manova, & Bisnovatyi-Kogan 1995, hereafter LRBK; Yi, Wheeler, & Vish- 
niac 1997) as well as the observed rotation-period distribution in YSOs (e.g., 
Konigl 1991; Edwards et al. 1993; Bouvier et al. 1993; CoUier Cameron 
& Campbell 1993; Yi 1994, 1995; Ghosh 1995; CoUier Cameron, Campbell, 
& Quaintrell 1995; Herbst et al. 2000). Furthermore, if the field is strong 
enough, then it may truncate the disk before it reaches the stellar surface 
and channel the inflowing matter to high stellar latitudes, where the (by now 
nearly free-falling) gas is stopped and thermalized in accretion shocks. This 
is the accepted explanation for X-ray pulsars (accreting magnetic neutron 
stars; e.g.. Lamb 1989) and DM Herculis stars (accreting magnetic white 
dwarfs; e.g., Patterson 1994), and it has also been proposed as the origin of 
the optical/UV "hot spots" in classical T Tauri stars (interpreted as accreting 
magnetic YSOs; e.g., Bertout, Basri, & Bouvier 1988; Konigl 1991; Edwards 
et al. 1994; Hartmann, Hewett, & Calvet 1994; Lamzin 1995; Bertout et al. 
1996; Johns-KruU & Basri 1997; Johns-KruU & Hatzes 1997; Martin 1997; 
MuzeroUe, Hartmann, & Calvet 1998). 

However, because of the complexity and diversity of the physical processes 
that are involved, the theory of the interaction between a magnetic star and a 
surrounding accretion disk remains incomplete. One key question is whether 
a steady-state description, which has been adopted in many of the models 
constructed to date, is appropriate. In a Keplerian disk around a star of mass 
M that rotates with angular velocity Q*, the gas interior to the corotation 
radius r^o = {GM /VtlY^^ rotates faster than the star, whereas the matter at 
larger radii rotates more slowly. If the star, the disk, and the magnetosphere 
that occupies the space between them can be treated as perfect conductors, 
then stellar magnetic field lines that thread the disk at any radius other 
than will undergo secular twisting. One possibility for establishing a 
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steady state is for the twisting of the field hnes to be countered by magnetic 
diffusivity in the disk (e.g., GL). However, in many cases a reahstic estimate 
of the disk diffusivity leads to values that are much too small to justify a 
steady-state description of a system with a dipolar field configuration except 
in the immediate vicinity of Vco (see § p.2| below). 

The behavior of twisted magnetic field lines that are anchored in a well- 
conducting medium has originally been considered in the context of the solar 
corona by means of semianalytic (e.g., Aly 1985; Low 1986; Wolfson 1994; Aly 
1995) and numerical (e.g., Mikic & Linker 1994; Amari et al. 1996a,b) tech- 
niques. However, direct insight into the field evolution in rotating accretion 
disks has by now been gained also from explicit investigations of magnetically 
linked star-disk systems, in which both semianalytic [e.g., van Ballegooijen 
1994 (hereafter VB); Lynden-Bell & Boily 1994; LRBK; Bardou & Heyvaerts 
1996 (hereafter BH)] and numerical (e.g., Hayashi, Shibata, & Matsumoto 
1996; Goodson, Winglee, & B5hm 1997; Miller & Stone 1997; Goodson, 
Bohm, & Winglee 1999; Goodson & Winglee 1999) approaches have again 
been employed. These studies have indicated that the applied twist leads to 
a strong expansion of the field lines away from the star. Initially, the evo- 
lution is quasi-static, with the azimuthal field component building up while 
the poloidal field structure changes relatively little, but past a certain point 
the expansion accelerates rapidly and in a finite time (corresponding to a 
total twist ~ vr) approaches a singular state characterized by the opening of 
at least a fraction of the field lines. 

The evolution of the magnetic field lines after the singular state is reached 
is another important question on which there has been no unanimity in the 
literature. According to one school of thought (e.g., LRBK), the opening 
effectively severs the magnetic link between the disk and the star once and 
for all. The opposite view is that the star-disk link is reestablished through 
rapid field-line reconnection in the magnetosphere, and that the field then 
continues to undergo successive episodes of twisting, expansion, and recon- 
nection (e.g., Aly & Kuijpers 1990, Goodson et at 1997, 1999). The latter 
picture is consistent with the results of numerical simulations of twisted coro- 
nal arcades, in which the addition of finite resistivity was found to result in 
reconnection across the current sheet that forms in that case (e.g., Mikic & 
Linker 1994; Amari et al. 1996a). It has, furthermore, been argued that, in a 
real magnetosphere, the current concentration will itself facilitate reconnec- 
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tion through the initiation of plasma microinstabihties that give rise to an 
anomalous resistivity (e.g., Lamb, Hamilton, & Miller 1993). To resolve this 
question in a quantitative manner, it is necessary to calculate the evolution 
of the magnetospheric current and particle densities as the field lines are be- 
ing twisted in order to determine whether an instability is likely to develop 
before the field expands to such a degree that it becomes effectively open. 

When the expanding, twisted field lines approach the singular state, they 
typically become sharply bent at the disk surface, giving rise to a strong radial 
stress that tends to push the field outward through the disk for reasonable 
values of the disk diffusivity. Previous treatments of the problem have either 
ignored this issue altogether (e.g., LRBK), concluded that it would lead to the 
effective expulsion of the field lines from the disk (BH), or else considered the 
intermediate case in which the radial field excursions average to zero (VB). 
The resolution of this issue depends in part on whether reconnection in the 
magnetosphere terminates the field-line expansion phase before the radial 
stress at the disk surface becomes very large, and it is thus directly related 
to the field-opening question. 

Even if the field configuration is not steady, the duration of the successive 
expansion/reconnection cycles is sufficiently short [of order the dynamical 
(rotation) time], that one could consider averaging the relevant quantities 
over the cycle period (cf. VB) and using them in modeling effectively steady- 
state star-disk systems (e.g., GL; Zylstra 1988; Daumerie 1996). However, 
the time-dependent nature of the magnetic interaction could be crucial to the 
origin of the observed accretion and outflow properties of disk-fed magnetic 
stars (e.g., Hartmann 1997; Goodson & Winglee 1999), and it may also help 
resolve some of the difficulties identified in steady-state magnetic accretion 
models (e.g., Safier 1998). 

In this paper we address various aspects of the basic issues listed above using 
a simple, largely semianalytic, approach. Our modeling framework, based on 
evolving a series of self-similar, force-free equilibrium magnetospheric field 
configurations, is described in § ^ In § ^ we clarify the condition for steady- 
state evolution in a disk with a given distribution of magnetic flux and dif- 
fusivity, and we then consider also the radial migration of flux across the 
disk. In § § we examine the opening of twisted magnetic field lines using 
a uniformly rotating, self-similar disk model, and we investigate the role 
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of reconnection and inertial effects in the magnetosphere above the disk in 
hmiting this process. In § ^ we generahze our analysis of the opening of 
magnetic field lines by using a non-self-similar model that incorporates a 
differentially rotating disk. In § ^ we discuss our results and relate them to 
recent numerical simulations. We summarize in § |^. 



2 Perfectly Conducting, Uniformly Rotating 
Disk Model 

In this section we describe a semianalytic model of a force-free magnetic 
field above a perfectly conducting thin disk. This model, which was first 
developed in 1994 by VBQ, is distinguished by its relative mathematical sim- 
plicity, and we use it to illustrate the relevant ideas and as a framework for 
our quantitative analysis. 

2.1 Self- Similar Configurations 

Following VB, we consider a uniformly rotating disk that is magnetically 
linked to a central, point-like star. This may be an adequate representation 
of the outer parts of a Keplerian disk, at radii r ^ r^o, where the beat 
frequency AQ [the difference between the rotation rate fid('") of the disk 
and the rotation rate fi* of the star] is almost independent of r (and is 
approximately equal to — We work in spherical coordinates {r,6,(j)) and 
assume that the system is axisymmetric. 

Since the disk is assumed to be infinitely conducting, the distribution of the 
poloidal magnetic flux \1/ on its surface is fixed and does not change as the 
star rotates relative to the disk: \l/d(r, t) = \l/d('"). This distribution serves 
boundary condition at the disk surface. 

^The self-similar model constructed independently, and at about the same time, by 
Lynden-Bell & Boily (1994), is mathematically identical to that of VB. This "spherical" 
model has been adapted to extended star-disk systems by BH, who, however, postulated 
different boundary conditions at the disk surface than VB. 
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We are looking for a self-similar solution, and, therefore, neither the bound- 
ary conditions nor the solution itself can possess a characteristic scale in r. 
Therefore, ^d(?^) must be a power-law function, 

*d(r) = -r-^ (2.1) 

Th 

where the coefficient C and the power exponent n are positive real constants. 

Our goal is to find a time-dependent solution for the magnetic field in the in- 
finitely conducting magnetosphere above the disk. Assuming that the plasma 
density is low enough for the Alfven speed va in the magnetosphere to be 
much greater than both the relevant Keplerian rotation speed and the 
sound speed Cg, the field configuration above the disk at any given time is 
given by a force-free equilibrium: 

V X B = a B, (2.2) 

where a{r,t) is a scalar function, constant along each field line. Since the 
boundary conditions of our 2-D problem [given by ^d('")] do not change with 
time, the only way time enters the problem is through the function a. 

Thus, the time-dependent solution is described by a sequence of equilibria 
parametrized hy t, or, equivalently, by the twist angle 

A$(t) = (f^d - ^*)t = Ant. (2.3) 



Following VB, the self-similar magnetic field can be written as 

C 



B = [Br, -Be, B^] 



„n+2 



sin^ 



(2.4) 



corresponding to 



\l/('r, 



C 



nr' 



9{0) 



where the functions f{0), g{0), and h{6) are also functions of time. Thus, the 
assumptions of axisymmetry and self-similarity enable the problem of finding 
3-D force-free equihbria to be reduced to a ID problem of determining the 
above three functions. 
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The boundary condition (2.1) implies 5'(7r/2) = 1, whereas the condition 
V • B = gives 

and the condition of self-similarity requires that a scale as 1/r, i.e., 

„(,,i) = ^!(M£M. (2.6) 



By integrating dr/dO = rBr/Be, the shape of the field line is 

r(^,*)=roW[^W]^/", (2.7) 

where ro(\I/) is the position of the footpoint of the field hne, labeled by the 
flux function ^, on the disk surface. 

Then, the condition that a is constant along the magnetic field immediately 
gives 

a(^,A$)=ao(A$)[^(^)]Vr (2.8) 
[Note that the function g{9) also depends on time, and therefore on A$.] 

Next, the 9 component of equation (2.2) enables h{9) to be expressed in 
terms of g{9): 

Finally, the (j) component of equation (2.2) gives the following second-order 
nonlinear differential equation for g{9): 

(^efl +n(n + lMe)^^^4mr^''--0. (2.10) 

The boundary conditions are g{0) — and g{'K/2) — 1. 

We note that equation (2.10) can also be obtained directly from the Grad- 
Shafranov equation for force-free magnetostatic equilibria: 

sm9 d ( 1 
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where the function F(^) in the nonhnear term on the right-hand side rep- 
resents the contribution of the toroidal field component and is related to 
q;(^, t) via 

a = F'{^). (2.12) 

This can be verified by making the self-similar ansatz (2.4) and using equa- 
tions (2.6)-(2.8). 

Equation (2.10) contains the single parameter oq, so one gets a one-parameter 
family of solutions that describe a sequence of equilibria. The dependence 
of g{9) on time is implicitly parametrized by ao(t). In order to find how ao 
changes with time, one needs to calculate the azimuthal twist angle A$ as 
a function of ao, and then use the relation (2.3). The twist angle A$ can 
be obtained by integrating the equation smOdcj)/ d9 — Bfj^/Bg along the field, 
yielding 

2.2 Magnetic Field Evolution 

For any given ao, one can solve equation (2.10) by direct integration and 
then use the obtained solution in equation (2.13) to derive the corresponding 
value of the twist angle. This is the procedure followed by VB. We have, 
however, found it advantageous to reformulate the problem so as to make 
the twist angle A$, rather than ao, the control parameter. In this way the 
time dependence of the sequence of equilibria becomes more transparent: 
for a given t, A$ is given by equation (2.3), and then the problem can be 
solved using A$ as the input parameter. This approach also has the merit 
of expediting the solution process. 

We have accomplished this goal by replacing the dependent variable g{6) by 
the function 4>{9), the twist angle of a given field line as a function of 9: 

figior-^.. (2.14) 

n + I Jo sm9 

The idea here is to express the function g{9) through (j){9) using equa- 
tion (2.14), and then substitute it into equation (2.10), thus obtaining a 
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differential equation for 0(^). Then the parameter A$ comes in through the 
boundary condition 0(7r/2) = A$, and, as we show, oq completely drops out. 



We have: 



-m] 



l/n 



1 



+ 1 '■■^ ' " sin 6 ' 
and after differentiating two more times one gets 



(2.15) 



k"2 



n 2 



ng sin 9 



ncf)' 



ng 



Using equation (2.10) to express g" in terms of g' and one finally gets a 
third order differential equation for 0(6'), which does not contain oq: 



6"^ 6' 

— (l-n) + — V(2cos2^-n)-0"(2n-l)/tan0-(n + l)0'^sin^0. (2.16) 
(p' sin 6 

The three boundary conditions are: 

1) m = 0; 

2) 0'(O) = — this can be used only if n < 2: indeed, near 9 = 0, g{9) 9'^, 
and so (j)'{9) ~ ^^/"/sin^ ~ ^(2/n)-i ^ as ^ if n < 2; 

3) 0(7r/2) = A$ — prescribed value. 

After the solution 0(6*) is found, one can use the boundary condition giji /2) = 
1 to determine Qq: 

ao(A$) = (n + l)0'(7r/2) . (2.17) 



In Figure |T] we plot ao(A$) for several values of the parameter n. As has 
been noted by VB, the dependence of on A$ is nonmonotonic, and for 
any value of between and ao,max = '^o(^'^'max) there exist two different 
solutions.0 

^Actually, there is an infinite series of bands of values of A$ where sohitions exist, 
separated by bands of forbidden values. Any two solutions with values of A(f> in a given 
band belong to the same topological class [e.g., the function g{9) has the same number of 
nodes, with the first band, < A$ < A$c, having zero nodes], whereas two solutions in 
different bands are topologically distinct and cannot be smoothly transformed into each 
other. 
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Table 1: PARAMETERS OF THE SELF-SIMILAR MODEL 



n 






'^0,max 


h 

'''d,ma,x 




1.0 


0.00 


1.23 


1.64 


0.82 


2.036 


0.5 


1.05 


1.41 


3.00 


2.00 


2.645 


0.25 


1.30 


1.48 


5.15 


4.12 


2.944 



A purely poloidal field (A$ = 0) is potential (oq = 0). As the field-line twist 
(A$) grows, increases and reaches a maximum at some finite twist angle 
A$jnax(^)- As A$ is increased even further, gq decreases and eventually 
vanishes at a certain critical twist angle A^df^) [> A$rnax(^)]- We refer to 
the part of the curve that corresponds to < A$ < A$rnax as the ascending 
branch, and to the portion corresponding to A$max < A$ < A$c as the 
descending branch. Both A$rnax and A$c are in general of the order of one 
radian and depend only on n (see Table |I]; A$o is the twist angle for which 
fin/2) = 0, see § pj). 

The nonexistence of a solution for ag > ao.max is in agreement with the result 
obtained by Aly (1984) for his Boundary Value Problem I, in which one 
prescribes values of the normal magnetic field component i?„ and a on the 
boundary of an infinite domain. We also draw attention to the change of 
curvature (from convex to concave) exhibited by the solution (for n < 1) 
near the critical point. The physical origin of this behavior is clarified in 
Appendix A. 

The behavior of the toroidal field component at the surface of the disk is 
given by the function hd = /i(7r/2, A$). According to equation (2.9), is 
equal to ao/(7^ + 1), so it also first increases with A$, reaches a maximum 
value /id,max = ao,max/(^+ 1) at A$ = A$max('T-) (scc Table 1), and then 
decreases to zero at A$ = A$c- The fact that twisting up the field lines does 
not result in an arbitrary increase of the azimuthal field component at the 
disk surface is central to our analysis of the evolution of resistive disks (§ 
This behavior is emblematic of twisted flux tubes in general. As noted in § |l|, 
such tubes initially exhibit a nearly linear growth of with A$, but after 
A$ exceeds ~ A^^ax the built-up magnetic stress causes the field to expand 
rapidly (in the meridional plane), with the field-line twist traveling out to 
the region of weakest field near the apex of the flux tube. The propagation of 
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Figure 1: The dependence of ao on the twist angle A$. 

the twist can be understood in terms of torque balance along the tube (e.g., 
Parker 1979). The evolution of the poloidal field components with increasing 
twist angle is illustrated in Figure 2. It is seen that, for A$ > A$max, the 
flux tube expands and becomes elongated, with the direction of elongation 
defining the apex angle 6'ap.Q 



2.3 Finite-Time Singularity 

As the critical twist angle A$c is approached (corresponding to — > 0), 
the solution of equations (2.10) and (2.16) blows up, with the field lines 
expanding to infinity and thus opening up (see Fig. |^). In this limit, the radial 
field component at the disk surface diverges even as the surface azimuthal 
field tends to zero. No equilibrium solution is found for twist angles > A$c- 
This behavior is generic to twisted flux tubes in 3-D and is characterized 
as a "finite time singularity" (e.g., Aly 1995): the magnetic field reaches a 

■^Formally, the apex angle 0ap is defined as the value of 9 where g{9) has a maximum. 
Since g{6) and r{9) are related via equation (2.7), this angle also corresponds to the point 
on the field line that is most distant fi'om the star. 
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n=0.5 



1 2 3 4 5 6 7 e 9 10 




1 2 3 4 5 6 7 8 fl 10 



3 4 5 6 7 8 9 10 



AO=0.0 AO=1.0 AO=2.5 

Figure 2: Meridional projection of the magnetic field lines for three values 
of the twist angle A$ in the case n = 0.5, plotted in the (arbitrarily chosen) 
interval re [1, 10]. 

singular state after being twisted for a finite time (or, equivalently, by a finite 
angle) . 

To analyze the asymptotic (oq — > 0) properties of the function g{9), we note 
that Oo can be rescaled out of equation (2.10) by the substitution: 

G{e) = g{e)a^, . (2.18) 

The equation for G{9) then becomes 

'""Tei^e'I) ^n(n+l)Gm + -^^lGm'-''" = 0. (2.19) 

The boundary conditions are: G{0) = and G{'k/2) = Oq. Thus, the pa- 
rameter Go has been transferred from the equation to a boundary condition. 
But now the transition to the limit ao — > can be easily made, since the 
solution of equation (2.19) does not blow up as the boundary condition at 
6 = n/2 approaches zero. We designate the solution of equation (2.19) with 
the boundary conditions 6*0(0) = G'o(7r/2) = as Go{0). This function 
depends only on n, and remains finite in the entire interval (0, 7r/2). The be- 
havior of the original function g{9) as the critical twist angle is approached 
is then obtained from 

g{e,ao,n) ^ Go{e,n)ao'' as ao ^ . (2.20) 
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Figure 3: Plots of the function Go{0) for three values of n. 

Figure ^ displays Gq{9) for several values of n. These solutions have a number 
of useful implications. 

First, by combining equations (2.11) and (2.20), one can calculate the critical 
twist angle A$c: 

7r/2 

A$c = ^ / gT{9)^ (2.21) 
n + 1 J sm U 



Second, it is seen from Figure ^ that the apex angle 6'ap, at which the function 
Gq{9) reaches its maximum, is very close to 60° for a wide range of values 
of n. In fact, 6'ap — ^ 60° as n — ^ 0, in agreement with the finding by Lynden- 
Bell & Boily (1994). The azimuthal current density in the magnet osphere, 
which is related to g{0) via equations (2.2), (2.6), (2.8), and (2.9), is also 
concentrated around this angle, and in the limit of very small n the current 
distribution collapses into a narrow layer that extends radially along 9 ~ 
(see Fig. As we discuss in § [4.3| , such a current layer is a natural potential 
site for rapid field reconnection in the magnetosphere. 

Third, similarly to the current density, the azimuthal flux, according to equa- 
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Figure 4: Plots of the azimuthal current density j(f){6) in the hmit A$ A$c 
for three values of n 



tion (2.9), becomes concentrated near the angle 6'ap. This concentration is 
associated with the outward propagation of the magnetic field twist, as we 



discussed in 5 2.2 



3 Steady-State Configurations 



In this section we consider the conditions under which a magnetically linked 
star-disk configuration can reach a steady state. We continue to assume, 
as in § ^, that the star as well as the magnetosphere can be approximated 
as perfect conductors, but we allow the disk to have a finite resistivity.^ 
This picture should be adequate for capturing the basic behavior of real 
systems. For ease of presentation, we use the self-similar model outlined 
in § D, which corresponds to a uniform value of the field-line twist angle 
A$ at any given instant. The self-similarity assumption in turn imposes 
a condition on the radial dependence of the vertically integrated electrical 



The role of resistive effects in the magnetosphere is explored in 
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conductivity E. However, since the conditions we derive are essentially local, 
we expect our basic conclusions to remain valid also in the more general case 
of a differentially rotating disk and a non-self-similar diffusivity. To further 
simplify the presentation, we first consider the case where the radial positions 
of the magnetic field lines in the disk remain fixed during the twisting process, 
so the diffusivity only affects the evolution of the azimuthal field component. 
We then examine the consequences of incorporating also radial field diffusion 
into the picture. 



3.1 Time Evolution of the Twist Angle 

We start by deriving a general expression for the time evolution of the twist 
angle in the presence of finite disk diffusivity. Focusing on a given disk 
radius r, we assume that the disk has a thickness 2H -C r, so that the 
normal field component at the disk surface is given approximately by i?d,z = 
Bz{z = H) ^ Bz{z = 0) = —Bg{r,TT/2) (where the approximate equality 
follows from the solenoidal condition on B). The azimuthal field component 
B^, which for a twisted field generally has a finite value B^^^ at z = H, is zero 
at 2; = on account of the reflection symmetry with respect to the midplane. 
Hence the disk carries a radial surface current density Kr ~ 2 H jr{z = 0), 
given by 

Kr = -cBd,4,/2TT . (3.1) 

Now, by Ohm's law, the radial electric field at the disk surface is 

Ed,r = Kr/T. - Vd,4,Bd,z/ c , (3.2) 

where Vd,^, = T^^d is the azimuthal speed of the disk matter and where S is 
related to the magnetic diffusivity t] through S ^ Hc^/2nr]. On the other 
hand, the azimuthal speed of the field lines that thread the disk at the given 
radius is fB,^ = — C-Ed,r/-Bd,z- The azimuthal speed of the field lines relative 
to the disk gas is, therefore, fB,(/) — Wd,^ = —cKr/TiB^^z ~ vBd,<i>/ HB^ ^- 

The time evolution of the twist angle is thus given by 

—— = An + ^ , 3.3) 

where, as in § |2|, AVt = fid ~ 
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The first term on the right-hand side of equation (3.3) represents the secular 
growth of A$ due to the differential rotation between the disk and the star, 
whereas the second term (in which c^/27rS f» ^/H) describes the azimuthal 
slippage of the field lines relative to the disk material brought about by 
the finite disk resistivity. One can alternatively obtain this equation by 
generalizing the circuit analysis previously used to derive the steady-state 
condition (see references below) to include an inductive contribution to the 
line integral in Faraday's law (i.e., c/E-dr = —d'^/dt). The result, in any 
case, is quite general and depends only on the disk being thin and on the 
magnetosphere being perfectly conducting; no other assumptions (such as 
mechanical equilibrium in the magnetosphere) need to be made. 

We expect that, in general, the ratio -Bd,<^/-Bd,^ in equation (3.3) will be a 
function of A$. In particular, if the magnetosphere is described by a force- 
free equilibrium, then this function is uniquely determined, and, for a given 
distribution of S(r), one obtains a closed equation for A$(r, t). For the 
purpose of illustration, we consider this equation in the context of the self- 
similar solutions discussed in § |^. In order for these solutions to be applicable, 
E(r) must scale as 1/r. In that case i?d,,^/-Bd,z(A$) = — /i(7r/2,A$) = 
— ao(A$)/(r;, + 1) [see eq. (2.9)], so equation (3.3) becomes 

^?^ = An-^^£i^. (3.4) 

dt 27rrE n+1 ^ ' 



By setting the time derivative in equation (3.3) equal to zero, one obtains the 
steady-state value of i?d,,^/-Bd,z, the azimuthal pitch at the disk surface. If the 
vertical field component in the disk is assumed to be fixed, this expression 
yields the azimuthal field component at the disk surface for which the twisted 
field configuration maintains a steady state, 

27rrSAr] 

-D</-,ss = -Dd,2 • (3.5) 



This result was previously obtained by Campbell (1992), LRBK, and BH 
using an electric circuit formulation. The same basic relation was used by 
GL to define the effective electrical conductivity of the star-disk system in 
terms of the time-averaged azimuthal pitch at the disk surface. In contrast 
to the GL approach, equation (3.5) focuses on the actual conductivity of the 
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disk and gives the unique value of the azimuthal pitch that corresponds to 
a genuine steady state. It is then possible to inquire whether any particu- 
lar system with a given distribution of disk surface conductivity S(r) and 
differential rotation rate AQ{r) will be able to attain a steady state. 

To answer this question, we take the magnetospheric field to be force free 
and consider the self-similar solutions described in § p.2[ Typically in these 
solutions, |i?d,,^(A$)| at first grows with increasing A$, reaches a maximum 
li?^^""! at A$max (with \B^^^''/Bd,z\ ~ 0(1) for n ~ 1) and then declines to 
zero at A$c- We now show that whether or not a steady state can be reached 
depends on the relative magnitude of -B™^'' and B^^ss- 

We define a maximum surface conductivity Smax by 



S 



max 



„2 Rmax 



27irAn Bd,z 





C2 




27rrAl] 



^d,max ) ("^-6) 



where the second equality gives the result for our self-similar model. If the ac- 
tual surface conductivity of the disk is large, S > Smax, then l-B™^""] < |-B</,,ss|, 
and there is no steady state. In this case the azimuthal resistive slippage is not 
strong enough to balance the twist amplification of the azimuthal field com- 
ponent, and a singularity is reached in a finite time (corresponding to A$c)- 
The effect of the resistive diffusion is merely to delay the onset of the sin- 
gularity but not to remove it. This is illustrated in Figure ^, which shows 
that the time delay At due to resistivity increases with decreasing E. The 
slowdown rate is largest near A$max because, for this twist, the azimuthal 
field |-Bd,(/>|, and hence the resistive slippage in the azimuthal direction, are 
maximized. As the critical twist A$c is approached, the time derivative of 
A$ reverts to its initial value [c?A$/c?t](A$c) = [dA^/dt]{0) = A^l because 
the azimuthal field at the disk surface vanishes at the singular point. 

Conversely, when the disk conductivity is small (S < Umax), -B^^"" > -B^^ss 
and the steady state described by equation (3.5) can be reached. If one starts 
with a purely poloidal potential field (A$ = 0, i?d,(/) = 0), then A$ initially 
grows linearly with time but subsequently levels off, tending to a steady-state 
value A$ss that is given by the condition 

5d,^(A<l>,,) =5,^,,,. (3.7) 

This situation is illustrated in Figure |[ 
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Figure 5: Time evolution of the twist angle for several values of the disk 
surface conductivity S above Smax (eq. [3.6]) in the n = 1 self-similar solution. 
The time t is given in units of 1/AQ and the surface conductivity in units 
of c^/27rAQr (Umax in these units is equal to /id,max = 0.82). 
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Figure 6: Time evolution of the twist angle in the n = 1 self-similar solution 
for the case in which the disk surface conductivity is less than Smax- The 
normalizations of t and of S are the same as in Fig. |^. 

Note that, although in principle two solutions are possible for a given value 
of l-Bd.^l < l-B™^""!, the steady-state solution described here always corre- 
sponds to the left (ascending) branch, i.e., A$ss < A$max (provided that 
one starts the evolution with the potential field A$ = 0). Also note that 
the steady state corresponding to the left branch is stable, whereas the one 
corresponding to the right (descending) branch is unstable in the following 
sense. Suppose we are at some steady state, and we decrease A$ slightly. 
Then, if the solution is on the left branch, |-Bd,<^| decreases and the azimuthal 
slippage speed becomes smaller than AQr, causing A$ to increase back to 
its original value and thereby restoring the initial state. Solutions on the 
right branch exhibit the opposite behavior: if A$ is decreased slightly, |-Bd,</)| 
will increase and the enhanced azimuthal resistive slippage will cause the 
twist angle A$ to drop even further, thereby moving away from the original 
state. 

We now address the question of whether real magnetically linked star-disk 
systems may be expected to reach a steady state. Astrophysical plasmas are 
typically highly conductive, so magnetic diffusivity is commonly ascribed to 
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an anomalous resistivity associated with fluid turbulence (e.g., Parker 1979). 
The turbulent diffusivity can be expressed as i] = /?t>turb-f^, where fturb is 
the speed of the largest turbulent eddies and /5 is a parameter < 1. The 
eddy speed could be determined by either thermal or magnetic processes; to 
maximize the value of rj, we set fturb = maxjcg , f a}- In magnetized accretion 
disks, the sound speed generally exceeds the Alfven speed (evaluated at the 
midplane) except well within the corotation radius. But this is also the 
region where the magnetic field rapidly becomes strong enough to enforce 
nearly complete corotation with the star, which, in turn, suppresses field-line 
twisting and inflation: this region is therefore not of much interest in our 
present discussion.^ Concentrating, therefore, on the case Cg 3> fA, we can 
write 



which is constant for an isothermal disk. Comparing equations (3.8) and 
(3.6), we obtain 

^^hrr • (3.9) 

This ratio is ^ 1 in a thin disk {H ^ r) threaded by a dipolar field (n = 1), 
except in the region immediately adjacent to Tcq. 

In the case of molecular disks around YSOs, it is possible to estimate r] 
directly from an explicit determination of the electron-molecule collision fre- 
quency for given temperature and density. Adopting the expressions given in 
Meyer & Meyer-Hofmeister (1999), we have r\ = lO^'^^T^^"^ {nn/ne) , where T3 
is the temperature in units of 10^ K, and where the electron-to-neutral num- 
ber density ratio is calculated by assuming ionization equilibrium of alkali 
metals (primarily potassium) and is given by log (ne/n„) = 6.48 — IO.94/T3 + 
0.75 log T3 — 0.51ogr;,„. As an illustration, we consider T Tauri stars, which 
are relatively slow rotators (mean rotation rate ^ 10~^ s~^), and for 
which we infer (assuming a 0.5 Mq star) r^o ~ 8.7 x 10^^ cm. D'Alessio 
et al. (1998) modeled accretion disks around such stars, and for a typical 
accretion rate of 10~^ Mq yr~^ and a disk "a parameter" of 0.01, we infer 

^We also note in this connection that the turbulence in the disk is often attributed 
to the action of the magnetic shearing instability, but that the latter ceases to operate 
(neglecting for the moment resistive effects) when ua comes to exceed ~ Cg in an unlinked 
disk, or ^ (csWd,^)^^^ in a magnetically linked system (e.g., Gammie & Balbus 1994). 
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from their results values of ~ 2 x 10^ K and ~ 7 x 10^^ cm~^ for the midplane 
temperature and particle density, respectively, at r^o- For these values, we 
get 7] ^ 7 X 10^° cm s~^, which is 5 orders of magnitude smaller than the 
nominal maximum turbulent diffusivity CsH. Although the estimated diffu- 
sivity increases at larger radii, we consider the region interior to r^o to be 
particularly relevant since the disk must extend to r < r^o if matter is to be 
accreted onto the star. Note also that, as r increases above Tco, \AVL\ — ^ 
and Smax itself decreases with r (oc 1/r). We conclude that the result (3.9) 
is likely to represent a lower bound on the ratio S/Smax in many practical 
applications. 

The preceding discussion indicates that it is unlikely that a disk with a dipole- 
like field configuration will achieve a steady state. It can, however, be seen 
from Table 1 that /id,max — ^-/n for n < 1. (Note in this connection that 
BH demonstrate that /ij— i>l/nasn— s>0.) Therefore, for sufficiently small 
values of n, Smax could be large enough for a steady state to be attainable 
even for realistic values of S. In § p.2| we argue that small effective values of n 
are also needed to achieve a steady state (at least in the time-averaged sense) 
if variations in the disk radial (and not only azimuthal) field component are 
taken into account. 

3.2 Radial Flux Diffusion 

We now consider the role of radial resistive diffusion. So far we have assumed 
that the radial distribution of flux in the disk does not change. However, 
since we include resistive diffusion in the azimuthal direction, we need for 
consistency to also consider diffusion in the radial direction and its effect 
on Sd,2(r). 

If the radial motion of the disk matter is vanishingly small, then the radial 
speed of the magnetic field footpoints in their resistive diffusion across the 
disk can be written as fB,r = cEd,(j)/Bd,z = cK^/T,Ba^z, where £^d,(/) is the 
azimuthal electric field at the disk surface and 2Hj^{z = 0) is the 

vertically integrated azimuthal current density (cf. eq. [3.2]). Setting = 
C-Bd,r/2vr, where i?d,r = Br{r,TT/2) is the radial magnetic field at the disk 
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surface (cf. eq. [3.1]), one gets 

VBAr, t) = /-^ = ^/(vr/2, ao{t)) . (3.10) 
/ttZj JDd,z 27rZj 

(It is interesting to note that in the case of an effective turbulent conductivity 
described by equation (3.8), we have VB,r{r,t) = /5cs/[vr/2, ao(t)], i.e., the 
footpoints move radially across the disk with the speed of order of the sound 
speed!) 



We have already noted in § |2]^ that -Bd,r/-Bd,z diverges as the twist angle 
approaches its critical value A<l>c- This can be formally demonstrated by 
noting that, in this limit, (see § ^3), and by using equations (2.5) 

and (2.20), which imply /(7r/2) = g'{7r/2)/n ao"Go(7r/2)/n as ao 0. 
But Gq{7i/2) is finite and independent of ao (see Fig. 3), so it follows that 

/(7r/2,ao)(xao", ao ^ . (3.11) 

Thus /[7r/2, ao(t)], and hence VB,T{i",t), diverge as A$ — >• A$c- 



Now, our results in § for the finite-time singularity were based on the 
assumption that the footpoints of the magnetic flux tubes are fixed in the disk 
throughout the field-line twisting process. In light of the fact that the radial 
diffusion speed diverges as the critical angle is approached, this assumption 
is not strictly valid, and in reality the footpoints will undergo some radial 
excursion Ar during the twisting time from A$ = to A$ = A$c- If Ar 
remains much smaller than r then our results will not change much. However, 
if Ar/r > 1, then our conclusions about the approach to a singularity will 
need to be reexamined. 

We can obtain a conservative estimate of Ar by working in the high-conductivity 
limit (S ^ c^/AQr, corresponding to a diffusivity rj <C AQHr). In this 
limit we neglect resistive slippage in the azimuthal direction and use equa- 
tion (2.3); this approximation becomes even better justified near the critical 
point, where Sd.t/- (see eq. [3.3]). Assuming that the surface conduc- 
tivity S does not change with time, the total radial footpoint displacement 
over the time interval [0, tc] is 
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^2 /-A^c 



/(7r/2,ao(A<l>))dA$ = 



2nEAQ Jo 

In the special case when E oc 1/r, resistive slippage does not break the self- 
similarity assumption. Indeed, in this case the radial displacement at any 
given radius scales as r, with the coefficient of proportionality increasing 
with twist angle but having the same value at all locations. Therefore, only 
the overall magnitude of the flux distribution is affected by the footpoint 
migration, but its self-similar scaling (and the value of the power-law index n) 
remain unchanged. 

We have already obtained the asymptotic behavior of /(7r/2) in the limit of 
small Oo as A$ — > A$c (see eq. [3.11]). The calculation of (dA^/dao) in 
this limit is more cumbersome and is given in Appendix A. On the basis of 
equation (All) from that Appendix, we can write 

A$ ~ A$c + ^Go . ao ^ , (3.13) 

which implies 

^<xar\ ao^O. (3.14) 
Using equations (3.11) and (3.14) in equation (3.12), we obtain 

/•A*c /■ao=0 / dA^\ 

J /(7r/2,ao(A$))dA$~ J ( __ j /(7r/2, ao)dao ~ 



/ 



ao=0 _ _ 

ag^ag (iao~logao, (3.15) 



which diverges logarithmically as ao — >^ 0. 

The result (3.15) demonstrates that, even if the conductivity is relatively 
large, one cannot neglect the radial field diffusivity in the disk as the critical 
point is approached. As A$ — > A$c; the radial displacement of the magnetic 
footpoints gets progressively larger and eventually reaches a level where the 
approximation 5^,2 (r) ~ const that was used in the derivation of equation 
(3.3) becomes inadequate. The outward motion of the footpoints would have 
the effect of decreasing Bd,r and preventing the surface radial field component 
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from blowing up. It can be argued, however, that this does not mean that 



the finite-time singularity discussed in § |2]^ would be avoided, although its 
onset would be delayed somewhat. The reason is that, as t — *• tc, only 
scales logarithmically with ao(t), but the radius r(ro,6'ap) of the apex point 
of a field line that remains anchored in the disk at some given value of ro 
increases as l/ao(t) (see eq. [4.28] below). [In that limit, ao{t) decreases 
with time as (tc — t)^''"; see eq. (All) in Appendix A.] This implies that the 
twisted field lines will expand much more rapidly in the magnetosphere than 
their footpoints will migrate inside the disk, so the conclusion that (barring 
inertial effects) they open in a finite time should continue to hold. Although 
inertial effects will intervene to reduce the field-line expansion speed in the 
magnetosphere to roughly the local Alfven speed (see § |C^ ), the latter will 
still be much higher than the field diffusion speed in the disk (which, for a 
turbulent diffusivity, is of the order of fturb = max{cs , ^a}; see § |0| ). The 
effect of radial field diffusivity will be even less of an issue if reconnect ion 
effects in the magnetosphere terminate the field-line expansion before the 



critical twist angle is attained (see § ^^), since the value of B^ j. in this case 
will remain bounded. 

In one proposed scenario (see § |I|), the twisted magnetic field lines reconnect 
as A$ approaches the critical value, and the process of twisting, expansion, 
and reconnection recurs in a periodic fashion. If this indeed is what happens, 
then it is interesting to inquire whether the system can maintain a time- 
averaged steady state despite the expected radial diffusion of the field lines. 
VB, who first addressed this question, proposed that this could happen if 
n < 1, since, in that case, the radial magnetic field component at the disk 
surface (and hence Ar) changes sign in the course of the field-line evolution 
(see Fig. H). Therefore, depending on the value of the twist angle at which 
reconnection occurs, there will be a value of n G (0, 1) for which B^^r averages 
to zero between the start of the twisting cycle (when the field is potential) and 
the instant of reconnection. An alternative possibility, which also requires 
n < 1, was proposed by BH for an exact steady state. It is based on the 
fact that a self-similar field configuration corresponding to < n < 1 attains 
Bd,r = at some twist angle A$o that is generally less than A^max (see Table 
1). The possibility then arises that A$o is equal to A$ss (which, according to 
the discussion in § ^A\ , is also generally less than A^max), so that a genuine 
steady state with B^,. = (and with the twisting due to the differential 
rotation balanced by the toroidal resistive slippage) is established. However, 
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for a given n, this is only possible for some special value of S, since A^^^ is 
directly related to S. For n < 1, the required value of S is unrealistically 
large (of order c^/rAQ; see § 

Both the VB and the BH proposals could probably be realized only in sys- 
tems in which n ^ 1. In the case of the VB picture, this follows from the 
expectation that the twist angle at which reconnection occurs is quite close to 
), whereas in the case of the BH suggestion it is a consequence 
of the fact that, for characteristic values of E, a steady state will typically 
not be established unless n is very small. We defer until §P a discussion of 
the plausibility of real astrophysical systems being represented by such low-n 
configurations. 



A$c (see § 4.3 



4 Effects of Plasma Inertia and Reconnection 
in the Magnetosphere 

As we noted in § |l], magnetic field reconnection across the current concentra- 
tion region that forms when a twisted flux tube begins to expand is a possible 
alternative to the opening of the field lines. Whether reconnection can com- 
pete with the opening process depends in large measure on the magnitude 
of the plasma resistivity in the region of current concentration. Classical 
collisional resistivity is generally too small to play a role, and one appeals 
to anomalous resistivity induced by current-driven microinstabilities. The 
instability criterion is that the current density exceed a certain threshold, 
which is typically proportional to the electron density n^. This suggests that 
the relevant quantity to examine is the ratio of the current density j^p to He 
(or p). Anomalous resistivity could thus be triggered not only by the current 
concentration near 6'ap but also by the drop in density brought about by the 
rapid expansion of the toroidal field as the critical twist angle is approached 
(see § pD.g 

In this section we explore this possibility by calculating the velocity and 

^It is very tempting to suggest that a similar mechanism (i.e., the triggering of anoma- 
lous resistivity through a rapid density decrease due to the expansion of magnetic field 
lines) may also work for coronal mass ejections in the sun. However, this does not seem 
to be likely, since the density in the solar corona is relatively large. 
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density evolution outside the disk using the sequence of force-free equihbria 
presented in § 0. As part of this analysis we also evaluate the inertial effects 
in the expanding magnetosphere, which must remain small for the force-free 
equilibrium approximation to be applicable. In order to isolate the resistive 
effects in the magnetosphere and to distinguish them from the corresponding 
processes in the disk (which were considered in § we assume (as in § ^ 
that the magnetic field remains frozen into the disk material at all times. 



4.1 Magnetospheric Velocity and Density Evolution 

We assume that the magnetosphere remains force-free throughout its evo- 
lution, with the plasma on the whole remaining a good conductor and re- 



sponding passively to the motion of the field lines. Subsequently, in § 
we address the question of when inertial effects become important. Then, in 
we address the issue of reconnection. 



We determine the rate of change of the density from the mass continuity 
relation, which requires us to first obtain the velocity field. Writing the total 
velocity as the sum of the parallel and the perpendicular components, 

V = V|| + v_|_ = f ||b -|- v_L (4.1) 

(where la is the unit vector along the magnetic field), we first calculate the 
perpendicular velocity v^. We do this by applying the ideal MHD flux- 
freezing condition to the given magnetic field configuration B(r,t). [We do 
not use the equation of motion for since the latter has already been 
employed (in lowest order, neglecting inertial and pressure terms) to obtain 
the underlying sequence of force-free equilibria.] 

Consider a surface of constant 6{= Oq) and a field line that at some moment of 
time t intersects this surface at a point with coordinates {r[9o,t],9o,(f)[9o,t]). 
After some infinitesimal time interval dt the same field line intersects the 
surface 9 = 9q at a. new point (r(6'o,t + dt),9o,(j){9o,t + dt)). We can thus 
define the vector u(^0;^) as representing the rate of motion of the point of 
intersection of the given field line with the surface 9 = const: 

u = (r(^,t),O,0rsin^) , (4.2) 
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where the dot symbol denotes an exphcit time derivative (d/dt). By the flux 
freezing condition, = u^. Using equation (2.7) to express r{9,t) in terms 
of g{9, t), one thus gets 

^-r = -[l ^)-r<^^^, (4.3) 

rgf sine -ghsiv? 9 

= — p2 ' (4-4) 

v^4. = — — + r0sm^ 1 — , (4.5) 



where 



/ I "R I 'n,+2 \ ^ 

P2(^) = / ' j sin^ Q^f sin^ Q + g^'iQ) + sin^ Q . (4.6) 




We now proceed to calculate the parallel component of the velocity, Vy . Un- 
like the perpendicular component, is determined from the equation of 
motion, whose component along b is given by 

where Fy is the acceleration due to the sum of the paraUel projections of 
all forces (including, in general, the centrifugal, Coriolis, gravitational, and 
pressure forces). Equation (4.7) implies that 

V\\ = -(v • ^)v\\ + Vx ■ ^ + . (4.8) 



We first calculate the parallel projections of the centrifugal and Coriolis 
forces. (These two inertial forces appear because we are working in a frame 
of reference that rotates with the angular velocity of the star.) We have 

i^centll = Feent " b = ^ siu Q\f {Q) siu^ Q + g{Q) COS Q\ (4.9) 

and 

i^corll = 2[v X • b = 2[vx X Jl,] • b = 
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-p- v±(i,{f sin^ 9 + g cos 9) — h sin 9{v±r sin 9 + v±0 cos 9) . (4.10) 

Next, we note that the gravitational force scales with distance as r~^, whereas 
the centrifugal and Coriolis forces are proportional to r. In this analysis we 
limit ourselves to large distances (r ^ Vco) in order for our self-similar model 
to be valid. Then wc can neglect the gravitational force. We also neglect any 
pressure forces. Thus we get 

i^ll =i^cent|| +i^cor||. (4.11) 



Now we are ready to calculate v\\. Using the fact that, due to self-similarity, 
I'll oc r, so that dv\\/dr — v\\/r, we can rewrite equation (4.8) as 

VrV\\ dv\\ dh /, X 

^11 = - + • 9^ + • • + ^11 ■ (4-12) 

We need to derive the expressions for the two terms involving the change in 
the unit vector b. First, using ■ b = 0, we rewrite • dh/dt in terms of 
the known functions f{9,t), g{9,t), and h{9,t) as 



dh 



(v±rf sin 9 + v±gg + v±^h sin 9^ . (4.13) 



Similarly, after some algebra one obtains 



. (v • V)b = [v^r{g"/n -g) + (v^e{n + 1) + v^^aog^/"") fsm9] + 

-p- [—v±rh sin 9 — v±gh cos 9 + v±^ (/ sin ^ + cot 9)] . (4.14) 

Equation (4.12), supplemented by equations (4.9)-(4.11), (4.13), and (4.14), 
is the partial differential equation that determines v\\(9,t). We solve it nu- 
merically by the finite differences method, assuming the following boundary 
conditions at 6' = and 9 = 7t/2: 

^(0,t)=0, ^(7r/2,t) = 0, (4.15) 
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and the initial condition at t = (corresponding to A$ = 0): 



v\\{e,0) = 0. (4.16) 

Once the velocity field is known for all times, the density evolution can be 
determined from the continuity equation, 

P = - V ■ (pv) = -l|-(rV,,) l-^^(smepve) . (4.17) 

or r sm 6 oo 

It is easy to see that, on account of the self-similarity of the velocity field, 
the evolution equation (4.17) preserves any power-law dependence of p on r. 
Thus, if we assume that p{r, t = 0) oc r"^, we can write 

p{r,e,t) =r-'Pp{e,t) . (4.18) 

The continuity equation (4.17) can then be expressed as 

V \ d 

p = -(3-p)p— — — {sinOpve). (4.19) 

r rsmd o6 

We solve this partial differential equation numerically with the initial condi- 
tion p{6, 0) = 1, and the boundary conditions 

|^(0,t) = 0, p{n/2,t) = l. (4.20) 

Figure ^ presents the results of our numerical computations for the case 
n = 1, p = 0, and Q^, = —1.0 (we choose the units of time so that AQ = 1; 
then our choice of Q^, corresponds to a star rotating with angular velocity 
equal to 1, while the disk is not rotating). It is seen that both the r and the 
6 components of the velocity vary smoothly as functions of 6, and that they 
grow rapidly as A$c is approached. The radial velocity component has a 
maximum at the apex point 6^^ ~ 0.98, whereas vq changes sign there.[] The 
plasma density around 6'ap decreases rapidly with 9, which corresponds to a 
significant amount of plasma being moved toward the symmetry axis, where 
it is concentrated in a narrow sector (A^ of a few degrees). As we discuss in 
§ ^, this mass concentration may be relevant to the formation of stellar jets. 

^Note that ^ap changes slowly with time, starting at 6'ap — tt/2 a,t t — and then 
gradually approaching its asymptotic value 9a.p{tc), which depends only weakly on n. 
Generally, 0ap(^c) is very close to one radian. Since in this section we are particularly 
interested in the behavior near t = tc, we consistently evaluate ^ap at that point. 
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Figure 7: Radial and 9 velocity components and the plasma density as func- 
tions of 6 at fixed radius at: (a) A$ = 0.5, 1.0, and 1.5, and (b) A$ = 1.5, 
2.0. These plots are for the n — 1 self-similar field configuration with a uni- 
form initial density distribution {p — 0). Time is measured in units of 1/AQ: 
in these units the critical twist angle A$c is attained at tc — 2.036. 
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Figure 8: The ratio of the current density to the plasma density p (in 
arbitrary units) as a function of 6 for A$ = 0.5, 1.0, and 1.5 in the n = 1 
self- similar model. 



Figure ^ shows the ratio of the current density to the mass density, which we 
consider to be a field reconnection diagnostic, as a function of 9. It is seen 
that, as the field lines expand, this ratio increases rapidly near 6 = ^ap, which 
suggests that the threshold for triggering anomalous resistivity (and, thereby, 
reconnection) could be exceeded. However, in order to reach a definitive 
conclusion, one needs to find out whether the threshold is reached before 
inertial effects become important. The answer to this question requires a 
more detailed knowledge of the behavior of the velocity and density at the 
apex in the limit t —>■ t^. We therefore perform an asymptotic analysis of the 
magnetospheric velocity field in this limit. 

The details of this analysis are given in Appendix B. Here we only summarize 
the main results by giving the asymptotic expressions for the velocity: 

v±r{r,9,p,t)c^-^ (4.21) 
vMr, ^ ^ap, t) ^ -r ^ + ^^~^7 , (4.22) 
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Table 2: ASYMPTOTIC VELOCITY AND DENSITY PARAMETERS AT 
THE APEX ANGLE 



n 


^ap Col^ap) 




Ai 


A2 


9i 


92 


1.0 


0.976 3.45 


1.99 


1.99 


11.96 


4-p 


-1-p 


0.5 


0.995 2.37 


3.87 


6.17 


36.4 


7.4-2p 


-0.4-2p 



and 



where 



(^1/" f) _ a 

v^^ir, 9 ^ t) ^ -r ' \y , (4.23) 



v^^{r,9,t)c^ri9-9,^)^, (4.24) 



A = ^ (3n + 4 ± V5n2 + 12n + s) (4.25) 



and K = Y 1 + G'o^"'(^^ap) /(^ + 1)^- For example, for n = 1 we have = k ~ 
2 and = 6k ^ 12. 

We can now proceed to calculate the behavior of the density at 6'ap as t ^ tc- 
Substituting the expressions for vj_ and v\\ into equation (4.19), we readily 
obtain 

pocr-P{t,-t)\ (4.26) 

where 

, ^ 3 — p n + 1 A 
q{n) = ^ + — . 4.27 

n UK, 

Table || lists the values of 6'ap, G'o(6'ap), n, A, and q for two representative 
values of n. Note that because we have a quadratic equation for A, two roots 
{Ai and A2) are possible for each n, giving rise to two different values (gi and 
^2) of the density power exponent q. It is not a priori clear how to choose 
the physically relevant value for A, and hence for q. However, based on the 
good agreement between the value of qi and the result of the full numerical 
solution for the case n = 1 (see Fig. 9), we conjecture that the appropriate 
root is generally the lower one, Ai. 

Figure ^ details the comparison between the predictions of the asymptotic 
analysis and the results of the direct solution of equation (4.19) for the 
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Figure 9: Log-log plot of p{9ap{tc),t) near t = tc for the n = 1, p = self- 
similar solution. The density is plotted in arbitrary units and the time in 
units of 1/ AQ. The dashed line represents the asymptotic solution p{9ap, t) oc 
{tc — t)^ derived in the text. 

case p = 0. It is seen that a reasonably good agreement is achieved when t 
approaches tc to within ~ 0.1/ AQ or so. Furthermore, by examining the plot 
in Figure |l^ of the entire time evolution of p(^ap(^c), i), one sees that there is 
a gradual transition in the interval t ~ 1.6 — 2.0 Af2~^ from the early linear 
dependence on t to the behavior [oc {tc — t^] predicted by equation (4.26). 
We note that a high value of q (as chosen in this figure) results in /)(^ap) 
becoming quite small well before the singularity. It is also interesting to note 
that the drop in density near t = tc is, in fact, faster than that corresponding 
to uniform expansion [guniform = (3 — p)/n]. This is due to the 6'-divergence 
of the flow near 6 = 6'ap. 

If, instead of fixing r, a given field line \i/ is considered, one can use equa- 
tions (2.7) and (B4) to obtain 

r(*,^ap,t) = ro(^)ao"^Gy"(M «^o(^)Ar]-^/"(tc-t)-'/", (4.28) 
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Figure 10: Normalized plasma density at the final apex point ^ap(^c) as a 
function of time (in units of 1/AQ) for n = 1, p = 0. 

and therefore, by equation (4.26), 

p(*, ^ap, t) oc po(*) Af)''+^/"(tc - oc {t, - t)i+4/n-A/. ^ ^4 29) 

where po is some typical initial density (at A$ = 0) on this field line. 

Having found the asymptotic behavior of the plasma density, we can now 
address the issues of inertial effects and reconnection in the magnetosphere. 
We start with the role of inertial effects. 

4.2 Inertial Effects in the Magnetosphere 

As the field lines expand near the singular point and the plasma velocity at 
the apex becomes larger and larger, there is a concern that it will become 
large enough to be comparable with the local Alfven speed, thus invalidating 
the equilibrium assumption. 

In order to investigate this question, we first need to find the behavior of the 
Alfven speed on an expanding field line at ^ = ^ap ast ^ tc- Consider a field 
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line \I/ intersecting the disk at the radius ro(\l/). Initially, at zero twist angle, 
the Alfven speed is large. 



C 

VA,o = VA{ro, ^ap, A$ = 0) ~ ^ AQ vq , (4.30) 

'^o vPo 



The inequality (4.30) justifies our underlying assumption that the magnetic 
field in the magnetosphere is in equilibrium. In fact, the ratio AQrQ/vAfi 1 
plays the role of a small parameter in the asymptotic analysis of the inertial 
effects in the magnetosphere. 

We now use equations (B7) and (4.28)-(4.29) to find the behavior of the 
Alfven speed near the critical point, 

VA{^,e.p,t) ~ t;^,o[Afi(te-t)]'/"-'^/'-^/'" ~ VA,o[mtc-t)f^'^-'/'. (4.31) 

(Interestingly, for = 1 one gets A = k, so va ^ const ~ va,o in this case.) 

Then, using equations (4.21) and (4.28), 



fr(6'ap,t)| roAQ 
VA{0^p,t) Va,0 



[Afi(tc-t)]"^^+"+^^ (4.32) 



It is seen that in this case the power-law exponent is independent of p and is 
always negative (for example it is equal to -2 for n = 1). This means that, 
sooner or later, the inertial effects on a given field line will become important. 
We can even estimate the time when this happens as 

2n 

An{t, - t)i„ ~ . (4.33) 



Note that inertial effects may, in fact, become important even earlier, but 
not near 6'ap. Indeed, on a significant part of a field line, p does not decrease 
as much as near 6'ap, whereas the typical velocities there are comparable 
to Vr(6'ap)- However, we shall use the estimate (4.33) anyway, because we are 
mostly interested in the neighborhood of 6'ap. 

The onset of inertial effects could in principle remove the finite-time singular- 



ity induced by the twisting of the field lines (see § 2^). However, we expect 
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that in practical applications the effective opening of the field lines will only 
be delayed by these effects, rather than eliminated altogether (see § || for 
a fuller discussion of this point). Furthermore, we anticipate that, even in 
the absence of inertial effects, the field-line expansion in real systems will be 
terminated by reconnection — although, as we show in the next subsection, 
inertial effects could also delay the onset of the latter process. 

4.3 Reconnection in the Magnetosphere 

In this subsection we address the prospects for reconnection in the disk mag- 
netosphere. We first derive the relevant equations in the self-similar model 
framework, and then use them to draw some quantitative conclusions using 
representative parameters for YSO systems. 

Reconnection may occur when the critical twist angle is approached. It is 
important to recognize, however, that the twisted field configuration does 
not give rise to a fully developed current sheet. Formally, a current sheet 
is characterized by an infinitesimal thickness in the limit 77 — 0, whereas 
in our self-similar solution the current is concentrated in a region of finite 
angular extent about 9 = 6'ap (see Fig. ^). The angular half-width of this 
region becomes very small only in the limit n — > 0, but for n ~ 1 it is of order 
~ 10°. Hence one can say that current is concentrated in a layer with an 
aspect ratio of order 10 or so; therefore, in order for reconnection to be of any 
significance, one needs a very low Lundquist number. In the framework of 
the Sweet-Parker (Sweet 1958, Parker 1963) reconnection model, one would 
need S = vp^L/r] ^ 100 (where L is the characteristic length of the current 
layer). As in almost any other astrophysical context, the required value of S 
is too small to be accounted for by the classical Spitzer resistivity. 

It can, however, be expected that near the critical point the ideal MHD ap- 
proximation will cease to be valid in the region of high current concentration. 
One interesting possibility is the development of an effective anomalous re- 
sistivity that is triggered by current- driven instabilities. For definiteness, we 
focus on the anomalous resistivity associated with the ion- acoustic /Buneman 
instabilities — the mechanism most commonly considered in the literature 
(e.g., Coroniti & Eviatar 1977; Galeev & Sagdeev 1984; Parker 1994). An- 
other possible mechanism for anomalous resistivity is lower-hybrid drift tur- 
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bulence initiated by a perpendicular current (see Krall & Liewer 1971 and 
Drake et al. 1984, and references therein). However, lower- hybrid insta- 
bility requires perpendicular electron drift, which is absent in the force-free 
configurations employed in our model. 

A current-driven anomalous resistivity is triggered when the current density 
exceeds a certain critical value, 



j = eriefd > jc = euev, 



c ) 



where Vd is the drift speed of the current-carrying electrons, and the criti- 
cal speed Vc is of the order of the thermal speed of either electrons or ions 
(depending on the nature of the instability) and is thus a function of the 
temperature T. The criterion for triggering reconnection can therefore be 
written as 

'->('-] iT) = '-^, (4.34) 
P \pJc Pemp 

where /ie is the molecular weight per electron in units of the proton mass. 

Assuming first, for simplicity, that the magnetosphere is isothermal (T = 
const), we infer from equation (4.34) that reconnection is most likely to start 
at the point where the ratio j/p is maximized. From our analysis we know 
that, at a fixed r, the current density is maximal, and the density is minimal, 
at the apex angle 9 = 9s.p (see Fig. ||). The radial position of the possible 
reconnection site cannot be determined within the framework of the self- 
similar model adopted here. In a realistic situation involving a Keplerian 
disk, however, we can argue that the ratio j/p will be greatest at the apex 
point of the field line that underwent the largest expansion, i.e., the field line 
with the largest twist angle. 

We now examine how the ratio j/p at the apex angle changes with time as 
t — > tc- According to equations (2.6), (2.8), and (B7), 

j{r,e^p,t) = -aB = --^Go (^ap)^^^— - c-^ ^^^^^ 

(4.35) 

where we reintroduced the speed of light c. If we fix the field line instead of 
the radius, we get 

rC 

Ji^^ ^ap, t) ~ ^:p3— [An{t^ - t)f- . (4.36) 



38 



Using equation (4.29) for p(\l/, 6'ap, t), we finally obtain the ratio j / p: 



P 



(^, ^ap, t) 



cC 



So Po, 



[A^{tc - t)] 



A/K-l-l/n 



(4.37) 



It is seen that the power-law exponent characterizing the asymptotic time 
evolution is independent of p, and, interestingly, is exactly equal to the ex- 
ponent in the asymptotic expression for vq that follows from equations (B9) 
and (B12). This implies that the time evolution of jif,/ p on a given field line 
is governed by the divergence of f ^ at = 6'ap. It is now also possible to 
estimate when the condition (4.34) for triggering anomalous resistivity will 
be satisfied, enabling reconnection to take place. The result is 



Af](tc - 1), 



nA/n—n—l 



PefTlr 



cC 



nA/K 



(4.38) 



Finally, by comparing the estimates (4.33) and (4.38), we can derive the 
condition for the onset of reconnection to occur before inertial effects become 
important in the magnetosphere. For example, for n = 1, we find 



Perrip cC V ^A,o 



< 1 . (4.39) 



So far in this subsection we have assumed that the temperature stays fixed 
during the evolution. However, because the field-line expansion is very rapid 
near the critical point, it might be more realistic to assume that the evolution 
of the plasma is adiabatic rather than isothermal. Indeed, let us consider the 
time scale for temperature equalization due to electron thermal conduction 
along the magnetic field lines. For the adopted fiducial parameter values (see 
below), the electron collisional mean free path is Ac — 10^ cm <C tq — 10^^ cm. 
The electron thermal conductivity is then Xe — ^cVth,e — lO^^cm^ s~^, and 
the temperature equalization time due to electron thermal diffusion can be 
estimated as ~ r^/xe — 10^ s ^ {AQ)~^. This motivates an examination 
of the adiabatic expansion limit .0 

^The thermal evolution of the magnetospheric plasma could also be affected by wave 
dissipation, heating by the stellar and disk radiation fields, recombination, etc., which 
complicate its precise determination. These effects may restrict the range of applicability 
of the adiabatic approximation to temperatures above some minimum value, which may 
be of order 10^ K. 
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Under the adiabatic assumption, the density drop near is accompanied 
by a drop in temperature T oc p^'^^^^ (where 7 is the adiabatic index), and 
by a corresponding decrease in the critical drift speed Vc oc oc ' " . 
One can thus write Vc ~ t'c,o(p/po)*''''~^^^^; where is the critical drift speed 
for the initial temperature Tq. 



The criterion for triggering anomalous resistivity can then be written as 



I 
P 



Po 



> 



(4.40) 



Correspondingly, we get 



AQ(ic ^)onset,ad 



(7+l)n(A/«-l)-47+2 



(4.41) 



By comparing the expressions (4.41) and (4.33), we derive the condition for 
reconnection to occur before inertial effects become important. For n = 1, it 
is given by 



/ic^p cC 

which generalizes equation (4.39). 



1 

27-1 



va,o 



<1, 



(4.42) 



To obtain quantitative estimates, we choose the following set of fiducial pa- 
rameters, which may apply to accreting, magnetized YSOs: 

= 1.5 X 10^^ cm; 
ro = 5i?* = 7.5 X 10^^ cm; 

B{R^, ^ = 0) = 1 kG =^ C = 1.7 X 10^6 Gcm^ (for n = 1); 
An = lO-^s-^ Uefi = 10W-=^; To = 10^ K; fx^ = 1. 

Using these values, we get: 

AQro = 7.5 x 10^ cm s"^; 

va,o = = 9 ■ 10^ cm 8"^ 

Afi(tc -t)in ~ 0.09 < 1; 

^c.o ~ Vth,e = 1-2 X 10^ cm s"^; 

Jo = cC/Anrl - 0.013 G s-^ 

^^d,o = Jo/ene,o ~ 27 cm s~^; 

An{t, - t)onset,iso ~ 3 X lO"" < An{t, - t)i„; 

AO(tc - ^)onset,ad ~ 4 X lO'^ < AVt{t^ - t),^ — for 7 = 5/3. 
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On the basis of these estimates we conclude that inertial effects should be- 
come important much earlier than the time when reconnection is first trig- 
gered. We note, however, that this conclusion depends on the rather uncer- 
tain value of the initial plasma density Uefi- if that density were lower, then 
reconnection could occur while the equilibrium assumption was still valid. 
For the isothermal case with n — 1, this requires (using eq. [4.39] and the 
adopted fiducial values of Tq, Tq, Bq and AQ) 

r^J^< 300cm-^ (4.43) 

corresponding to An(tc — t)onset,iso ^ 0.01, whereas for the adiabatic case 
(using eq. [4.42] and setting 7 = 5/3), the upper limit on the density is 

<J<lxl0^cm-^ (4.44) 

which imphes AD.{tc - t)onset,ad k, 0.03. 

Besides addressing the issue of the onset of anomalous resistivity, one has 
to consider whether the current-driven microinstabilities, once triggered, will 
in fact lead to an anomalous resistivity that is large enough to provide the 
required low value of the effective Lundquist number {S^s ^ 100). Using 
the Sagdeev (1967) estimate of the anomalous collision frequency associated 
with the ion-acoustic instability, 

u,ff = 0.01 Upiivjc,) {TJT,) (4.45) 

(where Cg is the speed of sound, uj-p is the plasma frequency, and the sub- 
scripts e and i refer to the electrons and ions, respectively), we obtain the 
corresponding anomalous resistivity from ?7cfr = c^mcZ/cfr/47mce^. Assuming 
that the electron and ion temperatures are equal, we then find the effective 
Lundquist number to be 

5e// = 2 X 10^ {uj,,L/c) {cjv^) , (4.46) 

where uj^ denotes the cyclotron frequency. We shall use — i'th,e even though 
the ion-acoustic instability can arise at lower values of the drift speed, since 
the effective resistivity is not significantly enhanced unless ~ Wth.c (e.g., 
Kadomtsev 1965). We then get an expression that, interestingly enough, is 
independent of both density and temperature, 

5e//^4x 103(u;,iL/c). (4.47) 
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The expression (4.47) is applicable once the condition (4.34) (or its adiabatic 
counterpart [4.40]) is satisfied. After the ion- acoustic microturbulence is 
initiated, the value of Ses at the apex of the expanding field lines decreases 
on account of the time evolution of a;ci oc i? ~ {C/tq) [AQ{tc — t)]^ and 
Lr^r^ ro/[An{tc - t)]. In fact, 

Seff ~ So An{t, -t)^0, t^t,, (4.48) 

where 5*0 = 4 x 10^u;ci{t — O)ro/c ~ 4 x 10^. This demonstrates that, even 
when excited, the ion-acoustic anomalous resistivity will saturate at a level 
that is too low to provide an efficient route for reconnection. 

It is, however, worth noting that, after being triggered, anomalous resistivity 
will be strongly localized. Over the past two decades, numerical simulations 
(Ugai & Tsuda 1977; Sato & Hayashi 1979; Scholer 1989; Erkaev et al. 2000) 
have shown that a strong local enhancement of resistivity may lead to a 
transition to Petschek's (1964) regime of fast reconnection. Nevertheless, 
it appears that, even if one assumes the maximum Petschek reconnection 
infiow velocity, Urec = va/^^S, the length of time available for reconnection 
(from the moment when anomalous resistivity is triggered until the singular 
point tc) will be too short for any significant amount of magnetic fiux to be 
reconnected. Indeed, the condition that all magnetic flux within the area 
of angular width ~ 0.1 around the apex is reconnected over the time 
interval {tc — t)rec can be written roughly as 

Urecito - t)rec > r = Tq ^^^^ ^_ ■ (4.49) 

Using Uj-ec = va/ InS* together with the expression (4.31) for va, we obtain 

[An{tc - t)rec]' = (^^\ [In So + Htc - t).ec] • (4.50) 

V ^A,0 J 

Neglecting the term ln{tc — t) compared with In 5*0 in equation (4.50), and 
taking A^ = 0.1, we get (tc — t)roc > 0.13/An, which greatly exceeds even 
our highest estimate [{tc — t)res,ad] for the available time. 

Our expressions for the onset time of ion-acoustic microturbulence (eqs. 
[4.38] and [4.41]) have been obtained on the assumption that the equilib- 
rium model remains applicable, which, as we have noted, requires the initial 
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densities to be sufficiently low {ricfl < 10^ — IC^cm"^; see eqs. [4.43] and 
[4.44]). We emphasize, however, that even if the initial densities are higher 
and inertial effects set in early enough to render our onset-time estimates 
inaccurate, we do not expect the evolution of the field lines to be affected in 
a qualitative way. Rather, we anticipate that inertial effects will merely delay 
the onset of reconnection. The actual triggering time of the microsinstability 
could well exceed the nominal critical time t^, but it may be expected to 
remain smaller than the effective opening time of the field lines, which will 
also be increased by the inertial effects (see § [4.2|) . In a similar vein, we 
expect the relative ordering of (t — tc)onset and (t — tc)rec (eq. [4.50]), and 
therefore our conclusions about the expected efficiency of reconnection, to 
remain unchanged by the time dilation induced by inertial effects. 

Finally, we note that a strongly developed MHD tearing turbulence could, in 
principle, provide an alternative route to reconnection (Strauss 1988). For 
the estimated aspect ratio of the current layer that arises in our model, the 
required fluctuation level that needs to be sustained in order for reconnection 
to be efficient is \6B/B\ ~ 0.1. It is unclear whether such a comparatively 
high level of turbulence saturation could be attained. Additional uncertainty 
is created by the expected suppression of the onset and growth of the tearing 
modes by line-tying effects (see Mok & van Hoven 1982); this issue has not yet 
been addressed in the literature and may provide an interesting direction of 
future research. Notwithstanding these caveats, we consider hyperresistivity 
produced by tearing-mode turbulence to be a promising mechanism for fast 
reconnection of the twisted field lines. 



5 Keplerian Disk 

The analysis presented in the previous sections employed the self-similar solu- 
tions presented in § ^. In this section we examine the more realistic case of a 
Keplerian disk and compare its behavior to that of the self-similar, uniform- 
rotation model. Since the Keplerian disk has a characteristic radial scale, 
namely, the corotation radius Vco, the self-similar semianalytic approach is 
clearly inapplicable. The problem becomes fully two-dimensional and re- 
quires numerical tools. 
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Figure 11: Geometry of the problem. 



To tackle this problem we developed a numerical code that enables us to find 
sequences of equilibria once the following two functions are specified on the 
disk surface. The first is the rotation law, i.e., the rate of relative twist as a 
function of radius along the disk surface, AQ{r). The second is the magnetic 
flux distribution on the disk surface, \l/d('")- Unlike in the self-similar case, 
these two functions do not have to be power laws. We first describe our 
numerical procedure and then present the results of our simulations. 

The computational domain consists of the outside of a sphere of radius i?* 
(see Fig. |TT|). In the remainder of this section, we normalize the radius r 
by i?*. Using the symmetry with respect the disk plane, we only consider 
the upper halfspace. 

To find the force-free magnetostatic equilibria, we solve the Grad-Shafranov 
equation 



The main difficulty in this equation is that the nonlinear term on the right- 
hand side [F F' is not given explicitly. Rather, it is determined im- 




(5.1) 
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plicitly by the rotation law via the condition 



(5.2) 



where 



1 de 



(5.3) 



Bor sin 6 



is an integral along the magnetic field line 

The time t in equation (5.1) is the parameter controlling the sequence of 
equilibria, and An(^d('")) is a prescribed function representing the rotation 
law. For example, for a Keplerian disk, one has = — (rco/T)^^^]. 

In the remainder of this section, we normalize t by 

Our goal is to find the time sequence of equilibria for a given rotation law. 
We start at t = with the potential field, which corresponds to F = 0. We 
then give a small increment to t (corresponding to an increment in A$ of the 
order of a small fraction of a radian), find the solution using the procedure 
described below, then proceed to the next moment of time t, etc. 

For each moment of time t we solve the system (5.1) -(5. 3) iteratively: at 
the k-th. iteration we use the result F('=)(^) of the previous iteration to ap- 
proximate the function [taking -F(^) for the initial guess {k = 0) to be the 
solution for the previous moment of time, or zero for t = 0], solve the ellip- 
tic equation (5.1), then use the solution ^'('^^(r, ^) to calculate the integral 
/{'=)( ^) along field lines, and then we update the function F{^) according to 



We repeat this procedure until the process converges. 

When solving the elliptic equation (5.1) for given t and k, we use the re- 
laxation method. We introduce a fictitious time variable r and then evolve 
^(r, r, 6) according to 



We first tried to solve this system using a uniform grid in spherical coor- 
dinates (r, 9) . However, in this case it is necessary to introduce an outer 



= Af](^)t//(^)(*). 



(5.4) 




(5.5) 
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boundary at some large radius r = rmax, where one encounters several seri- 
ous problems, such as the choice of boundary conditions and the treatment 
of the integral (5.3) for field lines that cross this boundary. To bypass these 
issues, we effectively place the outer boundary at infinity and use the trans- 
formation ^ 

x = ^, (5.6) 

which maps r = oo to a; = 0, while keeping the inner boundary (the surface 
of the star i?,, = 1) at x = 1. Correspondingly, we replace the uniform (r, ^) 
grid with a uniform (x, ff) one.^ 

After applying the transformation (5.6), equation (5.5) becomes 

This equation is integrated on a rectangular domain in the (x, 6) plane with x 
running from to 1 and 6 running from to 7r/2. There are four boundaries: 
the surface of the star x = 1, the axis 9 = 0, the outer boundary x = 0, and 
the surface of the disk 6 = 7i/2. On three of these the boundary conditions 
are particularly simple: 

^(x,^ = 0) = *(x = 0,^) = 0, (5.8) 

and 

^(x = 1,^) = ^,(^), (5.9) 

where the latter condition represents a prescribed magnetic flux distribution 
on the surface of the star, which is assumed to be infinitely conducting. The 
results we show correspond to a dipole field, 

^,{e) = smH, (5.10) 

normalized so that the total amount of flux through the stellar surface is 1. 

The boundary conditions on the fourth boundary (the equatorial plane 6 = 
7r/2) are somewhat more complicated because our model incorporates an 
inner gap between the disk and the star, which breaks this boundary into 

^We also tried the mapping x = 1/r but have found that x = works somewhat 

better because it allows one to pack more gridpoints at larger radii. 



46 



two pieces (see Fig. |TT|). Typically we place the inner edge of the disk at 
Tin = 1.5 (corresponding to = 2/3, \E'in = 2/3). The space inside the gap, 
1 < r < Tin, is filled with very tenuous plasma, just like the magnetosphere 
above the disk. Hence, the field lines crossing the equatorial plane inside the 
gap must be potential, and, because of the symmetry with respect to the 
midplane, the magnetic field inside r^^ has to be perpendicular to this plane, 

— {x>x;^,e = -) = 0. (5.11) 



In the region r > rin the magnetic field lines are frozen into the disk surface; 
thus the fiux distribution in this region is fixed, similar to the situation at 
the stellar surface: 

^{x<x;^,e = ^) = ^^{x), (5.12) 

where "^dix) is a prescribed function. In our illustrative examples we again 
employ a dipole representation, 

xij^(x) = x^ = - . (5.13) 



We start with the potential dipole field at t = 0. We then proceed through 
the sequence of equilibria by gradually increasing t (and therefore the twist 
angle), and using the solution for the previous value of t as the initial guess 
for the next value of t. Although our procedure can be used with any choice 
of Af2(r), our choice of the twist function was guided by the need for AQ{r) 
to vanish in the inner gap. For numerical convenience, we want this function 
to remain smooth as it approaches zero at Tin (which also makes physical 
sense, since we expect that near the inner gap the accreting gas will undergo 
a gradual transition from a Keplerian rotation law to corotation with the 
star). Along the rest of the disk surface, however, this function can be 
arbitrary. We investigated two particular cases: uniform rotation (Fig. |T^), 
wherein Af2(r) approaches a constant value as one moves away from the inner 
edge; and Keplerian rotation (Fig. [T^), where, for r > r^^, the rotation law 
becomes Keplerian with r^o = 6. 



We now turn to a description of our results. Figures |I3| and |1J present a series 
of contour plots of the magnetic fiux function for several instances of time 
for the uniformly rotating and Keplerian disk models, respectively. It is seen 
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that the basic behavior is very similar in both cases, the most important qual- 
itative feature being the rapid expansion of the field lines near ^ 60 — 70°. 
Just as in the self-similar model, this expansion is accompanied by a rapid 
rise of the toroidal current density j,^ in this region, as shown in Figures |1^ 



and 0. Figures and ^ describe the evolution of the function t) for 



the two models. 

In all the cases we studied, the evolution can be divided into two stages, 
distinguished by the time behavior of F(\E',t) on the field lines that have 
undergone the largest twist. For the Keplerian disk model, this field line is 
given by = 0.44, on which the absolute value of A$ is twice the asymptotic 
value at infinity (see Fig. |l2|b). Note that the function F(\l/ = t) serves as 
an indirect analog of ao(A$) in the self-similar model. The plot of F(\l/i,i:) 
is shown in Figure 

In the uniform-rotation disk model the analogy can be made more direct. A 
convenient way to describe the behavior in this case is to look at the evolution 
of the second derivative of Fi^^) at \I' = 0. Indeed, at large distances, r ^ rjn, 
the twist angle A$ approaches a constant (see Fig. [T^a), so one can expect a 
self-similar power-law asymptotic behavior of -F(^I/) in the limit \1/ ^ 0. Using 
a(\l/) = F'(\I^) and -F(O) = 0, one can express F(\l/,t) using the notation of 
nias 

F{^,t)=Jad^ = J^-^. (5.14) 





In the case n = 1 considered in our numerical calculations one has ro(\l/) = 
1/^^, so 

F(vl.,t) = ^^, v,^o. (5.15) 

This quadratic behavior is indeed exhibited by our calculated solution. 

The foregoing considerations have motivated us to select the time evolution 
of 

for a direct comparison with the corresponding function in the self-similar 
model. Figure ^ demonstrates that on the ascending branch of the solution 
the agreement is very good, but that on the descending branch there is a 
significant deviation that is reflected in different values of the critical twist 



«oW = 3-j;^k=o (5.16) 
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Figure 13: Sequence of magnetic field contours for the uniformly rotating 
disk model. 
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Figure 14: Sequence of magnetic field contours for the Keplerian disk model. 
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Figure 15: Toroidal current density as a function of ^ at a given radius 
(r — 10) at different moments of time for the uniformly rotating disk model. 
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Figure 16: Toroidal current density as a function of at a given radius 
(r = 10) at different moments of time for the Keplerian disk model. 
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Figure 17: for the uniformly rotating disk model. 
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Figure 18: Fi^^t) for the Keplerian disk model. 
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angle A$c- This discrepancy can be attributed to the fact that the innermost 
field lines have smaller twist angle and therefore are not infiated as much 
as in the self-similar case. Consequently, the magnetic stresses driving the 
expansion are somewhat weaker and the opening of the field lines is delayed. 
Still, we see from Figure ^ (as well as Fig. ^) that the basic behavior is 
the same. As t is increased, ao(^) (and F(\l/i,t) in the Keplerian case) at 
first rises until it reaches a maximum at some instant tmax, and subsequently 
starts to decrease, just as in the self-similar case. During the first stage 
(the ascending branch, t < tmax) the shape of the field lines does not change 
significantly, but during the second stage (the descending branch, t > tmax) 
there is a rapidly accelerating expansion of the field lines, which approach 



an open state (see Figs. |T3| and |T^. This quahtative behavior appears to be 
universal and is independent of the details, such as the particular rotation 
law of the disk. However, certain quantitative features, such as the values 
of ^max and ao,max [or -^(^'i, tmax)] do depend on the particular parameters 
(e.g., Tin, Tco, ^^1, etc.). 

The two evolutionary stages are also characterized by distinct convergence 
properties of the numerical solutions. During the first stage the convergence 
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t 

Figure 20: Evolution of the function ao(t) [eq. 5.16)] for the uniformly rotat- 
ing disk model. The solid line represents the result of the numerical calcu- 
lations discussed in this section, whereas the dashed line shows the behavior 
of the corresponding function in the self-similar model discussed in § ^ 
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is rapid and robust, but during the second stage it progressively slows down 
and one has to update the function F(^) more frequently, and finally we 
reach a point where we need to terminate the computation. The reason for 
this is that the field lines have by that time expanded very strongly, and our 
iteration process no longer converges at a given resolution. Instead, we get 
an unphysical reconnection in the middle of the domain, which renders the 
numerical procedure inconsistent. Although increasing the resolution helps, 
the computation time necessary for convergence grows dramatically.[^ For- 
tunately, by the time we are forced to stop the computation the field lines 



have already expanded by a large factor (see Figs. |2T| and ^), so an ex- 
trapolation of the functions ao(''^) and F(\E^i,t) becomes possible. We deduce 
that these functions reach zero at a finite twist angle (about 2.7 rad for the 
uniform-rotation case and 2.0 rad for the Keplerian case), which corresponds 
to opening of the field lines in a finite time, similar to the behavior of the 
self-similar solutions. 

We thus see that a good case can be made for the finite-time singularity oc- 
curring not only in the self-similar solution (see § |2.3| ) but also in a broader 
class of models. To strengthen this point, we present a simple physical argu- 
ment showing why this should be the case.Q 

Consider the small-F limit of equation (5.1), in which the nonlinear term 
FF'(\1') for some chosen field line is smaller than, say, the first linear term 
on the left-hand side at typical distances of order the footpoint radius ro(\E'). 
By a dimensional analysis, this implies 

(5.17) 

Small values of F may correspond to two qualitatively different solutions. 
The first solution is close to the potential field, with the nonlinear term being 
unimportant everywhere on the field line. The other solution corresponds to 
greatly expanded field lines. In the latter solution, a given field line stretches 
out to distances so large that the linear terms (dimensionally proportional to 

^"^This behavior can be understood from the fact that our relaxation procedure represents 
a diffusion process with a diffusion coefficient Dx oc x^, and the diffusion time over a 
distance fa ~ a; is proportional to /D^ ~ x^'^ — r^, which diverges as the field lines 
expand. 

^^For a different, more mathematical argument in favor of a finite time-singularity in 
sheared, axisymmetric, force-free magnetic fields, see, e.g., Aly 1995. 
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Figure 21: The expansion factor as a function of magnetic flux in the uni- 
formly rotating disk model. 
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Figure 22: The expansion factor as a function of magnetic flux in the Kep- 
ler ian disk model. 
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\l//r^) become small and the nonlinear term gives an important contribution 
to the equilibrium structure of the distant portion of the expanded field. For 
a given value of F, we can define a characteristic radial scale rp where this 
state of affairs will prevail, 

M^,'^) = ;p^»^oW. (5.18) 

The expression (5.18) also gives an order-of-magnitude estimate of the posi- 
tion rap(\I/,t) of the apex of the field line \1/ at the time when -F(\l/) has the 
value F(^,t). 

Our next step is to use this estimate to calculate the twist angle corresponding 
to the field line \1/ at time t. This can be easily accomplished by using 
equations (5.2) and (5.3) and assuming that there are no thin scales (i.e., 
current sheets) in the 6 direction. This means that the part of the field line 
that corresponds to r ~ rj? is finite in 6. Then one can write: 

/HO 1 
^ Bgr{^,e} sm {6} {Berj^nm 



Now, (i?er)min = —[{d'^/dr)/ sin^^]mm can be estimated simply as (\I'/r) 
^'/rap(^',t) ~ ^/r^. Thus, 



A$(^, t) ~ F(^, t)^ ~ 0(1). (5.20) 



Therefore, to lowest order in F, the twist angle becomes independent of F 
as F ^ for a given field line. That is, as F goes to zero and the field lines 
open, A$ approaches a finite value and we have a finite-time singularity. 



6 Discussion 



In this section we consider the implications of our work to the various mod- 
eling issues discussed in § |l|. We start with one of our key results, namely, 
the finding (§ |^) that real astrophysical systems are unlikely to reach an 
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exact steady state unless they have significant negative radial field compo- 
nents at the disk surface in their untwisted (potential) state (corresponding 
to an n ^ 1 magnetic field configuration in the self-similar representation). 
We now argue that such low-effective-n configurations are not likely to be 
realized in axisymmetric, magnetically linked star-disk systems. 

Although models in which the stellar magnetic field is pulled into the disk 
have been considered in the literature (e.g., Ostriker & Shu 1995), and i?d,r 
is < even for a strictly dipolar field if the disk has a finite height, it 
is difficult to imagine that strong negative radial field components would 
naturally develop in the course of disk accretion onto a magnetized star. In 
fact, simple global disk models (e.g.. Spruit & Taam 1990) indicate that, so 
long as the disk is at least partially conductive, it will tend to exclude the 
(dipolar) stellar magnetic field and cause some of the magnetic flux to remain 
trapped inside its inner radius, thereby giving rise to a positive radial field 
component at the disk surface. Furthermore, there has been no indication 
from any of the numerical simulations carried out to date that the magnetic 
field in a disk interacting with a central dipole can evolve into an n <^ 1 
configuration. Although the published simulations have all been strictly 2- 
D, such systems are even less likely in 3-D. This is because the steady-state 
values of |i?d,(^i/-Bd,x| are expected to be well in excess of 1 for n <^ 1 (see Table 
1), implying that the field configuration might be susceptible to a disruptive, 
nonaxisymmetric instability (e.g., Kadomtsev 1966). 

It is conceivable, however, that negative values of B^ r would occur in the 
case of a nonaxisymmetric field configuration. Nonaxisymmetry could be in- 
troduced by the misalignment of a stellar dipolar field axis with the disk axis, 
which is, in fact, the common interpretation of many of the observed period- 
icities in accreting magnetic stars (e.g., Lamb 1989; Patterson 1994; Bertout 
et al. 1988). Alternatively, it could arise from the presence of isolated stellar 
magnetic loops that thread the disk (e.g., Bertout et al. 1996; Safier 1998). 
The latter possibility is consistent with the inferred (by radio observations) 
existence of extended (> 10 i?*), organized magnetic structures in certain 
YSOs (e.g., Andre et al. 1992). The axisymmetric self-similar model with 
n < 1 could provide a qualitative representation of the field evolution in this 
case, but some important aspects, such as the possibility that only a fraction 
of the field lines open first (e.g., Amari et al. 1996b), or the development 
of kink and/or buckling instabihties (e.g., Parker 1979), might be captured 
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only in a fully 3-D calculation. In any case, although such systems might 
attain a steady time-averaged configuration, they could never be in a strictly 
steady state. 

Another new result of our work has been the demonstration (§ p.2|) that, 
even in disks with a comparatively low diffusivity, in which the field lines are 
twisted until they open up, one cannot consistently neglect the radial migra- 
tion of the magnetic footpoints in the disk as the critical twist angle A$c is 
approached. We have argued that it is unlikely that a state with B^ j. = 
(even just in a time-averaged sense) can be attained unless, again, the mag- 
netic field configuration corresponds to that of an n <^ 1 self-similar model. 
Under these circumstances, a disk that undergoes cycles of successive opening 
and closing (by reconnection) of twisted field lines (and, to a lesser extent, 
even a disk in which a steady-state twist is maintained) should manifest a 
secular evolution of the radial flux distribution. We have exphcitly shown, 
however, that the rate of outward flux diffusion in the disk is signiflcantly 
smaller than the rate of expansion of the twisted fleld lines in the overly- 
ing magnetosphere. This provides a strong argument against the suggestion 
made by BH that the fleld lines in low-diffusivity disks remain perpendicular 
to the disk surface at all times. BH's conclusion that the twisting process will 
tend to expel to the outer edge of the disk almost all the stellar fleld lines 
that initially threaded it can be questioned on similar grounds. Although 
an outward flux migration might take place, we consider it unlikely that it 
would lead to a flux rearrangement in the disk that is as drastic as envisioned 
by BH. In fact, in certain steady time-averaged flux distribution may 

be established. This is evidently the case in the time-dependent model pre- 
sented by Goodson et al. (1999) and Goodson & Winglee (1999), in which 
inward radial mass motions that counter the outward radial flux diffusion 
take place in the disk. 

We have also investigated (§ §) the density redistribution in the magneto- 
sphere brought about by the motion of the twisted fleld lines. We showed 
that, as the critical twist angle is approached, the density in the apex region 
of the elongated flux ropes drops precipitously, with some of the depleted 
mass being pushed toward the rotation axis, where it becomes concentrated 
in a narrow angular sector about ^ = 0. A similar axial density enhance- 
ment has been found in the numerical simulations of Goodson et al. (1997, 
1999; see also Hayashi et al. 1996), where it is directly associated with the 
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formation of an overdense, "knotty" jet. As already noted by Goodson et 
al. (1999), this mechanism of producing axially outflowing condensations 
requires a good alignment between the rotation and magnetic dipole axes. 
Although our sequence-of-equilibria representation cannot capture the com- 
plex dynamical evolution of the simulated jets, it is gratifying that this simple 
model nonetheless allows a quantitative description of the basic process that 
generates their density structure.0 

Although inertial effects slow down the field-line expansion, one can expect 
that the twisted field lines will still effectively open in a finite time. By this 
we mean that for any region of a given size threaded by a finite amount of 
flux, the apexes of the majority of the field lines will leave the region un- 
der condsideration after being twisted by a finite angle. This is indeed the 
behavior found in numerical simulations. For example, in the model of Good- 
son et al. (1999), even though inertial effects are enhanced by a relatively 
high initial mass loading of the magnetosphere and by subsequent episodes of 
post-reconnective mass ejection from the disk, the field lines typically open 
in less than one differential star-disk rotation period (as stated in Goodson 
& Winglee 1999), which is to be compared with the time (< 1/2 period) pre- 
dicted by our model in the absence of inertial effects (see § ^ and § The 
basic question is then whether the field lines eventually reconnect, resulting 
in a repetitive cycle of infiation and reconnection, or whether the opening 
process is essentially a one-time affair. 

The 2-D numerical simulations carried out to date do not provide a clear an- 
swer to this question. Although reconnection was found to take place in the 
work of Hayashi et al. (1996), this was probably a consequence of their adop- 
tion of an unrealistically low threshold value of j/p for the onset of anomalous 
resistivity and of a rather strong functional dependence [t] oc (j/p)^] above 
the threshold. Numerical diffusivity was evidently also responsible for the 
reconnection events observed in the simulations by Romanova et al. (1998) 
of twisted, disk- anchored magnetic loops. In the case of the numerical model 
of Goodson et al. (1999), the nominal value of the numerical resistivity- 

^^The simulations by Hayashi et al. (1996) and Goodson et al. (1997, 1999) also show 
an equatorial density enhancement induced by mass being pushed away from the apex 
regions of the expanding flux loops toward lower latitudes. This feature of the density 
evolution is not reproduced by our calculations on account of the fixed-density boundary 
condition that wc imposed at the disk surface (see eq. [4.20]). 
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based Lundquist number in the innermost zone was ~ 100 at the start of 
their run. As noted in § [4.3| , our self-similar model results indicate that a 
value in this range is compatible with current-driven Sweet-Parker reconnec- 
tion. However, the reconnection that actually took place in this simulation 
was triggered by inward mass motions at the inner radius of the disk that 
strongly pinched the magnetosphere (see also Miller & Stone 1997). 

Our analysis of the 2-D self-similar model has led us to the conclusion that 
ion-acoustic microinstability will be triggered too late to result in significant 
reconnection before the field lines effectively open. The crux of the prob- 
lem is that the current layer that forms in this geometry along the field-line 
elongation direction is comparatively wide (aspect ratio ~ 10 for an initially 
dipolar field), so the effective Lundquist number that is required for effi- 
cient reconnection is rather low and is attained only relatively late in the 
field-line opening process. Tearing-mode turbulence could, in principle, lead 
to reconnection on the required scale if the tearing fluctuation level is suffi- 
ciently high. As discussed by Strauss (1988), such turbulence can occur in 
3-D, giving rise to hyperresistivity that may result in fast reconnection. We 
can, in fact, argue more generally that, in contrast with the idealized 2-D 
configurations that we employed in our model, real 3-D systems may exhibit 
nonaxisymmetric behavior that could significantly broaden the range of pos- 
sible scenarios. For example, it seems quite probable that, as the field lines 
are twisted, the 2-D equilibrium will become unstable to nonaxisymmetric 
MHD perturbations, such as the kink mode (a possibility already mentioned 
by LRBK). The detailed MHD stability analysis that is required for a closer 
examination of this suggestion is beyond the scope of this paper, but it seems 
likely that such instabilities could induce a transition to a new, nonaxisym- 
metric configuration that might have real (i.e., infinitesimally thin) current 
sheets. Such current sheets could form sites of rapid reconnection that would 
enable the system to relax to a less stressed state, and it is quite probable 
that the field-line expansion would then terminate much sooner than in the 
2-D scenario considered above. This possibility could be checked by means 
of future 3-D numerical simulations. 

We thus regard it as highly plausible that reconnection, leading to periodic 
(or quasi-periodic) evolution, would be the realistic outcome of the field-line 
twisting process. This inherently nonsteady behavior enhances the attrac- 
tiveness of this scenario as a modeling framework. In particular, in the case 
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of YSOs, it could be relevant to the interpretation of such phenomena as 
X-ray flares, the formation of a "knotty" axial jet and of a less collimated 
and slower (but more massive) disk outflow, and the heating of the permitted 
line-emitting gas in the magnetosphere (e.g., Hayashi et al. 1996; Hirose et 
al. 1997; Miller & Stone 1997; Hartmann 1997; Shu et al. 1997; Goodson et 
al. 1997, 1999; Goodson & Winglee 1999). Furthermore, since the successive 
openings and reconnections of the field lines have been inferred to modulate 
the gas inflow rate onto the star and the net torque exerted on the star by 
the accreting matter, they would also have a direct impact on the mass and 
spin evolution of the protostar. It may nevertheless be possible to construct 
stationary, global models of accretion disks in such systems by obtaining a 
suitable time average of the magnetic field configuration. This is, in fact, the 
approach originally undertaken by GL and subsequently elaborated upon by 
Zylstra (1988) and Daumerie (1996). In particular, Daumerie developed a 
model in which the flux distribution along the disk is calculated by approxi- 
mating the magnetic field in the magnetosphere as being current-free and by 
using an effective conductivity (based on the value of B^^^/ B^^^ through a 
relation analogous to our eq. [3.5]) to evaluate the poloidal field components 
in the disk. We note, however, that inasmuch as the azimuthal currents 
(which determine the disk contribution to the poloidal magnetic field) flow 
entirely inside the disk, one should employ the actual disk conductivity in 
this calculation instead of the above ansatz, as the latter could significantly 
overestimate the correct value in regions where the effective surface conduc- 
tivity exceeds the limit (3.6). It may also be possible to refine the treatment 
of the magnetospheric field by using the approach outlined in § |^ to calcu- 
late the magnetic field configuration in the fully force- free formulation (i.e., 
without neglecting current flows in the magnetosphere). After identifying 
the time at which reconnection is most likely to occur, one could proceed to 
derive the time-averaged values of B^^^ij") / B^^zij) and BA,,f,{r) / B^^^zij") and 
use them as boundary conditions for the disk model. 

Past numerical simulations have proved invaluable in clarifying some appar- 
ent difficulties and in pointing the way to further theoretical progress. For 
example, Safier (1998) has raised concerns about the ability of the magnetic 
field lines to thread the disk rather than be swept out in a thermally driven 
stellar wind. However, the numerical simulations (which typically do include 
thermal effects) have demonstrated that a field with an initially dipolar topol- 
ogy can maintain a configuration in which it continuously threads the disk 
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even as some field lines open and drive an outflow, and that various time- 
dependent effects (such as field-line twisting and mass motions in and out of 
the disk plane) combine to increase its strength above naive, steady-state es- 
timates. Nevertheless, more refined simulations (affording better control on 
the field diffusivity and incorporating 3-D effects) are needed to fully resolve 
this and other still-open issues. 



7 Summary 

Disk accretion onto a magnetic star occurs in a variety of astrophysical con- 
texts, ranging from T Tauri stars to X-ray pulsars. The magnet ohydro dy- 
namic interaction between the stellar field and the accreting matter can have 
a strong effect on the disk structure, the transfer of mass and angular momen- 
tum between the disk and the star, and the production of bipolar outffows. 
In this paper we concentrated on a key element of this interaction — the time 
evolution of the magnetic field lines that thread the disk resulting from the 
relative rotation between the disk and the star — and studied it using sim- 
plified models of axisymmetric, force-free fields. In particular, we employed 
the semianalytic similarity solution first derived by van Ballegooijen (1994) 
and Lynden-Bell & Boily (1994) to construct a sequence of magnetospheric 
equilibria parametrized by the relative twist angle A<I>, but we tested the gen- 
erality of the basic conclusions of this model (which strictly applies only to 
a uniform relative angular speed) by constructing numerical solutions of the 
Grad-Shafranov equation for systems with rotation laws that approximate a 
Keplerian disk. Our main results can be summarized as follows. 

• On the assumption that both the star and the magnetosphere can be 
approximated as being perfectly conducting, the behavior of the twisted 
field lines depends on the magnitude of the surface conductivity E of 
the disk. A steady-state configuration can be established only if, at a 
radius r in the disk, S(r) < Sinax(^) ~ c^/r|Ar2(r)| (eq. [3.6]). Weakly 
ionized protostellar disks as well as disks characterized by a turbulent 
diffusivity typically do not satisfy this condition except very close to the 
corotation radius, implying that the field lines undergo secular infiation 
and will open up when a critical twist angle A$c (~ 2 rad for an initially 
dipolar field configuration) is reached. 
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If inertial effects in tlie magnetosphere could be neglected, then the 
field lines would open in a fraction of the differential star-disk rotation 
period. However, inertial effects intervene, slowing down the expansion 
speed at the apex of the twisted field line to the local Alfven speed 
(which, for an initial dipolar field and a power-law density distribution, 
remains constant throughout the evolution in the self-similar model). 
This will delay the effective opening of the field lines, but will not by 
itself prevent it from occurring. 

If E(r) > Tiuiax{r) and the field lines undergo a strong expansion, then 
the radial magnetic stresses at the disk surface will cause the field 
lines to start migrating outward. A similar conclusion was reached by 
Bardou & Heyvaerts (1996), who went on to suggest that the field lines 
will be rapidly expelled from the disk. We argued, however, that this 
is unlikely to happen in view of the fact that the radial diffusion in the 
disk is much slower than the field-line expansion in the magnetosphere. 
If the field lines expand and reconnect in a periodic manner then a 
time-averaged steady state in which the mean surface radial field is 
zero is in principle possible (van Ballegooijen 1994), but this situation 
will not arise if the initial field configuration is dipolar (although it 
might conceivably pertain to isolated stellar fiux tubes that thread the 
disk). 

The opening of the field lines induces a redistribution of the magneto- 
spheric density p. In particular, it leads to a strong depletion of mass 
near the apex of the expanding field lines (which typically elongate in 
a direction making an angle > 60° to the rotation axis) and to a pro- 
nounced density enhancement near the rotation axis. The decrease in 
density along the direction of elongation (which also marks the location 
of a current layer) is conducive to the triggering of microinstabilities 
(such as ion-acoustic), since the instability onset condition corresponds 
to a threshold value of j / p. The density enhancement near the axis 
(which is largely a byproduct of the axisymmetry assumption) repre- 
sents a possible mechanism for the formation of an axially outflowing 
condensation (generalizing to a "knotty" jet if the process is repeti- 
tive), as suggested previously (e.g., Goodson et al. 1999) on the basis 
of numerical simulations. 

The question of whether the evolution of the twisted field culminates 
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in the effective opening of the field hnes (e.g., Lovelace et al. 1995), or 
whether it is a repetitive process (e.g., Aly & Kuijpers 1990) of field- 
line expansion and reconnection (not necessarily back to exactly the 
same field configuration that existed before the most recent expansion, 
in which case it will not be strictly periodic), is still open. The re- 
connection events observed in the numerical simulations carried out to 
date are by and large due to numerical diffusivity and therefore do not 
provide conclusive answers. Based on an analysis of the 2-D self-similar 
model, we concluded that anomalous resistivity derived from the ion- 
acoustic/Buneman instability will be triggered too late for significant 
reconnection to take place across the comparatively wide (aspect ratio 
~ 10 for an initially dipolar field) current layer that forms in this case 
before the field lines effectively open. We suggested, however, that hy- 
perresistivity associated with tearing-mode turbulence (Strauss 1988) 
might lead to efficient reconnection, although whether such turbulence 
could arise in the configuration under consideration is not clear. In ad- 
dition, it is conceivable that even faster reconnection could take place 
as a result of 3-D kinking of the twisted field lines that would create 
genuine (i.e., infinitesimally thin) current sheets. This issue might be 
resolved by MHD stability calculations and 3-D numerical simulations. 
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A CALCULATION OF {dA<!> /dao) NEAR THE 
CRITICAL POINT 

In this Appendix we present a more detailed analysis of the behavior of the 
magnetic field near the critical point. Our main goal is to find the asymptotic 
behavior of (dA^/dao) near A$c- We use the results of this analysis in § |3]^ 
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to determine the effects of resistive diffusion near the critical point, and 
also in § ^ in studying the role of inertial effects and reconnection in the 
magnetosphere. 

First, we note that calculating ((iA$/(iao) is not a straightforward task, 
because it requires the knowledge of the deviation of g{9) from the asymptotic 
solution G'o(6')ao"'- Indeed, if one simply uses the lowest order expression 
g{9) ~ G'o(6')ao", then, upon substituting it into equation (2.13), one gets 
equation (2.21). We see that in order to find the behavior of (dA^/dao) near 
the critical point, one needs to know the next-order (in a^) correction to g{6) 
as ao — > 0. 

Let us write 

g(9) = G(5)ao " = Go{9)aQ'' + SG{e)aQ'' , (Al) 

where G{9) is the solution of equation (2.19) with boundary conditions 
G(0) = 0, G{7i/2) = < 1. 

We want to find the correction 6G{6). In order to do this, let us introduce a 
small parameter e. We choose e so that 5G{ti /2 — e) ^ Gq{ti /2 — e) and at 
the same time e <^1. As can be seen from Figure ^, the derivative dGo/d9 
at = 7r/2 is in general a finite (meaning that it does not scale with ao) 
number that depends only on n. Thus, Go{ti /2 — e) ~ edGo/d9{n/2). Next, 
SG{7i/2 — e) is of order SG{7i/2) = ciq, so we must choose e in the range 

a[} < £ < 1 . {A2) 

Inside the e- vicinity of = 7r/2 we can approximate sin^ ~ 1 = const, so we 
get 

G"i9) + n{n + 1)G + ^^G^+'Z" = . (A3) 

with G{n/2) = a[J. 

In this region G'(^^) ^ 1, so we can neglect the nonlinear term compared with 
the other terms, and hence the solution can be written as 

G'inner(^) ^ ^0 + (^-^/2) -const. (A4) 

In the region 9 < 7t/2 ~ e, the correction 5G is small compared with Go{9), 
and we can linearize equation (2.19). We then obtain the following linear 
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equation for 5G: 

'"^K +»(» + l)^Gm + ^*GW[G.«))lV» = 0. (^6) 

Actually, near 9 — t:/2 the last term of equation (A5) becomes negligible 
compared with the other terms (because Gq{9) <^ 1 there). Since this is the 
only term that distinguishes this equation from the one describing the inner 
region 9 > 7r/2 — we can ignore the difference and extend the range of 
applicability of equation (A5) all the way up to 6* = tt/2, where we set the 
boundary condition SG{7t/2) = a^. 

Thus, the solution of the linear equation (A5) with the boundary conditions 

5G{0) = 0, 5G{7r/2) = {A6) 
provides accurate (to the lowest order in oq) correction to the function g{9). 

Using the function SG{9) we can now determine the asymptotic behavior of 
A$(ao) in the hmit ao ^ 0. We rite 

A$ = A$c + SA^ , {A7) 

where 

SA^ = / S[G'/^{9)]-^ . {A8) 

In order to estimate this integral, let us choose some ei ~ ao~^ <^ 1, < 
X < n. Then, within the ei-vicinity of ^ = n/2 we have G{9) ~ Oq + Gq{9) 
and Go{9) ~ [dGo/d9]\g=T,/2iO — 7r/2) ~ Ciei. We can thus estimate that in 
this region |5[G^/"(^^)]| < 6*2^}^", where Ci and C2 are positive constants of 
order 1. 

Correspondingly, we estimate the contribution to the integral (A8) from this 
region as 

< CseJ+'Z" , {A9a) 
where C3 is a positive constant of order 1. 



— / S[G'/-{9)] 

+ 1 Jir/2-ei 
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In the rest of the integration domain, we can write 5[G'^/"(^)] ~ (l/n)(5G)Gl/"' 
which yields 

1 r^/2-ei , AO 1 /■7r/2-ei 1 , ^0 

n + lJo ^ ^ '^sm9 n + 1 Jo n " ^ ' smO ^ ' 

Since bG ~ Oq, the contribution from this region is of order oH^. We see that, 
under the additional condition that 

e^+V" « «n ^ ^ < ^ ^ (^10) 

this contribution is much larger than the contribution from the ei-vicinity of 

Q — which we therefore can neglect. Thus, wc immediately sec that the 
correction to the twist angle A$ is of the order a^. Now let us determine 
the coefficient. 

Notice that, as one takes the upper limit of the integral (A9b) to 6 — 7r/2, 
i.e., as ei 0, the integral converges as a function of ei. This is be- 
cause, as ^ ^ 7r/2, 5G ^ a"^ = const and G'y""^(^) ~ (7r/2 - 
so r'''-''5GGy''-\e)de/sme ~ r/'-^H7r/2 - Of'^-^de ~ e}/" ^ as 
ei 0. 



Therefore, to lowest order in ao, we can write 

^ ^— - / 5GGt-\e)^ = e(n)a^ . (.Ill) 

n(n + 1) JO sm U 

Using the functions 5G and Gq{6) obtained above, we get, for example, 

(5A<l>(n = 1, oo ^ 0) ft! -0.17 ao 

and 

(5A$(n = 0.5, ao ^ 0) = -0.22 a^^ . 



B ASYMPTOTIC ANALYSIS OF THE VE- 
LOCITY FIELD NEAR THE APEX AN- 
GLE 

In this Appendix we present an asymptotic analysis of the velocity field near 
the apex angle ^ap in the limit t ^ te- 
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First we do some preparatory work. Near the critical point t = the behavior 
of g{9) in the vicinity of the apex angle ^ap can be described by 

g{9, ao) ~ ao-"Go(^) = "Go(^ap) [1 - l^^^^^y^] , (Bl) 

where we used equation (2.20) and defined n = n{n) = —Gq{9i^p) / Go{9g,p) . 
Using equation (2.19), 

/. = n(n + l) + (S2) 

Next, in order to get the time dependence, we need to relate ao to time. 
Using equation (All) we can write: 

""^^^sL^^ih:- ^^^^ 



Therefore, 



Furthermore, on account of equations (2.5) and (2.9), 

n sm n sm AO f — i^. 



Zap #1. o±±± i/ap i-ij 1^ f t-c 

and 



(S5) 



, ao"G'o"^^'^"(6'ap) ^ G'o"^^^"(6'ap) C 1 /o^x 

^^P^ (n + l)sin^ap (n + l)sin^apAQi-ie ■ ^ ^ 

In addition, similarly to A$ = AO, we can write (j){9ap,t) — AO^/^, where 
( is some constant of order 1. 

Finally, we derive the expression for P{9s_p,t) (eq. [4.6]): 

P(^ap, t) = Kg{9,p) = /,Go(^ap)^^ , {Bl) 



where « = ^1 + G'o/"(^ap)/(n + 1)^ . 
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Substituting these expressions into equations (4.3)-(4.5) for the components 
of v_L, we get 

v±r[y^p) - r - -- - — - , {B8) 

(-0, 



and 



We can now calculate the parallel velocity using the equation of motion (4.12). 
The centrifugal and Coriolis terms in this equation can be neglected in the 
vicinity of the critical point t = tc because their contributions are, respec- 
tively, of zeroth and first order in {t — tc)^^, whereas the contributions of 
most of the other terms are of order [t — tc)~^ (and therefore dominate in 
the limit t ^ tc). 

Substituting the expressions (B4)-(B10) for g, f, h, P, and Vj_ into equa- 
tions (4.12)-(4.14), we obtain an equation for v^i{9 ~ 6'ap,t): 

. _ n + 1 e - e^p dvii f||n + 3 9 ~ O^p {n + l){n + 2) 

^" ~ ~l^~de ^ ~^ t-t^ 86 ^ T^t-U ~ ^'^{t-t^y ■ 

(511) 

Equation (Bll) implies that v\\ can be written as 

v\\{r,9,t)=r{9-9,p)j^, (512) 

where A — A{n) is a numerical constant. After substituting this expression 
into equation (Bll), we get a quadratic equation for A: 

A' - Ak{3 + 4/n) + ^2 (^ + l)(^ + 2) ^ ^ _ ^^^3^ 

Equation (B13) has two solutions for each value of n: 

^1 2 = — f 3n + 4 ± V5n2 + 12n + s) . (514) 
2n ^ ' 

For example, for n = 1, we find A\ — k'^2 and = 6/« ~ 12 . 
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